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Abstract 

We study the evolution of long-period comets by numerical integration of their orbits, a more realistic 
approach than the Monte Carlo and analytic methods previously used to study this problem. We follow the 
comets from their origin in the Oort cloud until their final escape or destruction, in a model solar system 
consisting of the Sun, the four giant planets and the Galactic tide. We also examine the effects of non- 
gravitational forces and the gravitational forces from a hypothetical solar companion or circumsolar disk. 
We confirm the conclusion of Oort and other investigators that the observed distribution of long-period 
comet orbits does not match the expected steady-state distribution unless there is fading or some similar 
process that depletes the population of older comets. We investigate several simple fading laws. We can 
match the observed orbit distribution if the fraction of comets remaining observable after m apparitions 
is oc m~ o e 1 (close to the fading law originally proposed by Whipple 1962); or if approximately 95% of 
comets live for only a few (~ 6) returns and the remainder last indefinitely. Our results also yield statistics 
such as the expected perihelion distribution, distribution of aphelion directions, frequency of encounters with 
the giant planets and the rate of production of Halley-type comets. 

1 Introduction 

Comets can be classified on the basis of their orbital period P into long-period (LP) comets with P > 200 yr, 
and short-period (SP) comets with P < 200 yr; short-period comets are further subdivided into Halley-type 
comets with 20 yr < P < 200 yr and Jupiter-family comets with P < 20 yr (Carusi and Valsecchi 1992). 
The boundary between SP and LP comets corresponds to a semimajor axis a — (200) 2 / 3 AU = 34.2 AU, 
which is useful because: (i) it distinguishes between comets whose aphelia lie within or close to the planetary 
system, and those that venture beyond; (ii) an orbital period of 200 yr corresponds roughly to the length of 
time over which routine telescopic observations have been taken — the sample of comets with longer periods 
is much less complete; (iii) the planetary perturbations suffered by comets with periods longer than 200 yr 
are uncorrelated on successive perihelion passages (see footnote 2 below). 

LP comets are believed to come from the Oort cloud (Oort 1950), a roughly spherical distribution of 
some 10 12 comets with semimajor axes between 10 3 5 and 10 4 5 AU. The Oort cloud is probably formed from 
planetesimals ejected from the outer Solar System by planetary perturbations. LP comets — and perhaps 
some or all Halley-family comets — are Oort-cloud comets that have evolved into more tightly bound orbits 
under the influence of planetary and other perturbations. Jupiter-family comets probably come from a quite 
different source, the Kuiper belt outside Neptune, and will not be discussed in this paper. 

The observed distribution of the ~ 700 known LP comets is determined mainly by celestial mechanics, 
although physical evolution of the comets (e.g. fading or disruption during perihelion passage near the Sun) 
and observational selection effects (comets with large perihelion distance are undetectable) also play major 
roles. The aim of this paper is to construct models of the orbital evolution of LP comets and to compare 
these models to the observed distribution of orbital elements. 

This problem was first examined by Oort (1950), who focused on the distribution of energy or inverse 
semimajor axis. He found that he could match the observed energy distribution satisfactorily, with two 
caveats: (i) he had to assume an ad hoc disruption probability k = 0.014 per perihelion passage; (ii) five 
times too many comets were present in a spike (the "Oort spike" ) near zero energy. Since most of the 
comets in the Oort spike are on their first passage close to the Sun, he argued that they may contain volatile 
ices {e.g. CO, CO2) that create a large bright coma for the new comet, but are substantially or completely 
depleted in the process. When the comet subsequently returns (assuming it has avoided ejection and other 
loss mechanisms) , it will be much fainter and may escape detection. Most of the decrease in brightness would 
occur during the first perihelion passage, and the brightness would level off as the most volatile components of 
the comet's inventory are lost. This "fading hypothesis" has played a central role in all subsequent attempts 
to compare the observed and predicted energy distributions of LP comets. 

In § ^ we examine the observed distribution of LP comet orbits. The basic theoretical model of LP comet 
evolution is reviewed in § 0. The simulation algorithm is described in § m and the results are presented in 
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§ ||. The simulations and results are described in more detail by Wiegert (1996). 

2 Observations 

The 1993 edition of the Catalogue of Cometary Orbits (Marsden and Williams 1993) lists 1392 apparitions 
of 855 individual comets, of which 681 are LP comets. The catalog includes, where possible, the comet's 
osculating orbital elements at or near perihelion. When studying LP comets it is often simpler to work with 
the elements of the orbit on which the comet approached the planetary system (the "original" elements). 
These can be calculated from the orbit determined near perihelion by integrating the comet's trajectory 
backwards until it is outside the planetary system. Original elements are normally quoted in the frame of 
the Solar System barycenter. Marsden and Williams list 289 LP comets that have been observed well enough 
(quality classes 1 and 2) that reliable original elements can be computed. 

The difference between the original elements and the elements near perihelion is generally small for the 
inclination perihelion distance q, argument of perihelion u> and longitude of the ascending node f2; when 
examining the distribution of these elements we will therefore work with the entire sample (N — 681) of LP 
comets. The original semimajor axis and eccentricity are generally quite different from the values of these 
elements near perihelion, so when we examine these elements we shall use only the smaller sample (N = 289) 
for which original elements are available. 

2.1 Semimajor axis 

The energy per unit mass of a small body orbiting a point mass Mq is — ^GMq/o, where a is the semimajor 
axis. This is not precisely the energy per unit mass of a comet orbiting the Sun — the expression neglects 
contributions from the planets and the Galaxy — but provides a useful measure of a comet's binding energy. 
For simplicity, we often use the inverse semimajor axis x = 1/a as a measure of orbital energy. The boundary 
between SP and LP comets is at x = (200 yr)~ 2 / 3 = 0.029 ALT 1 . 

Figure [j] displays histograms of x = 1/a for the 289 LP comets with known original orbits, at two different 
magnifications. The error bars on this and all other histograms are ±1 standard deviation (a) assuming 
Poisson statistics (a = iV 1 / 2 ), unless stated otherwise. 

The sharp spike in the distribution for x < 10~ 4 AU -1 (the "Oort spike") was interpreted by Oort (1950) 
as evidence for a population of comets orbiting the Sun at large (a > 10 000 AU) distances, a population 
which has come to be known as the Oort cloud. Comets in the spike are mostly "new" comets, on their first 
passage into the inner planetary system from the Oort cloud. 

2.2 Perihelion distance 

Figure || shows the number of known LP comets versus perihelion distance q. The peak near 1 AU is due 
to observational bias: comets appear brighter when nearer both the Sun and the Earth. The intrinsic distri- 
bution N(q) (defined so that N(q)dq is the number of detected and undetected LP comets with perihelion 
in the interval [q, q + dq\) is difficult to determine. Everhart (1967b) concluded that N(q) oc 0.4 + 0.6q for 
q < 1 AU, and that for q > 1 AU, N(q) is poorly constrained, probably lying between a flat profile and one 
increasing linearly with q. Krcsak and Pittich (1978) also found the intrinsic distribution of q to be largely 
indeterminate at q > 1 AU, but preferred a model in which N(q) oc q 1 / 2 over the range < q < 4 AU. 
Shoemaker and Wolfe (1982) estimated N(q)dq cx 500g - 175 for q > 1.3 AU. 

These analyses also yield the completeness of the observed sample as a function of q. Everhart estimates 
that only 20% of observable comets with q < 4AU are detected; the corresponding fraction in Shoemaker 
and Wolfe is 28%. Kresak and Pittich estimate that 60% of comets with q < 1 AU are detected, dropping 
to only 2% at q = 4 AU. Clearly the sample of LP comets is seriously incomplete beyond q — 1 AU, and 

2 Angular elements without subscripts are measured relative to the ecliptic. We shall also use elements measured relative to 
the Galactic plane, which we denote by a tilde i.e. i, Q and u. 
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Figure 1: Distribution of original inverse semimajor axes of 289 LP comets at two different magnifications 
(panels a,b) and for the 170 LP comets with the most accurate (Class 1) orbits (panels c,d). Data taken 
from Marsdcn and Williams (1993). There is no obvious difference between the top and bottom panels, 
suggesting that observational errors in the inverse semimajor axes are unimportant. 
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the incompleteness is strongly dependent on q. In comparing the data to our simulations we must therefore 
impose a g-dependent selection function on our simulated LP comets. We shall generally do this in the 
crudest possible way, by declaring that our simulated comets are "visible" if and only if q < q v , where q v is 
taken to be 3 AU. This choice is unrealistically large — probably q v — 1.5 AU would be better — but we find 
no evidence that other orbital elements are correlated with perihelion distance in the simulations, and the 
larger cutoff improves our statistics. We shall use the term "apparition" to denote a perihelion passage with 
q < q v - 

We have also explored a more elaborate model for selection effects based on work by Everhart (1967a,b; 
see Wiegert 1996 for details). In this model the probability p v that an apparition is visible is given by 

[ if q > 2.5 AU, 

py(q) = I 2.5 - (q/1 AU) if 1.5 < q < 2.5 AU (1) 
[ 1 if q < 1.5 AU 

The use of this visibility probability in our simulations makes very little change in the distribution of orbital 
elements (except, of course, for the distribution of perihelion distance). For the sake of brevity we shall 
mostly discuss simulations using the simpler visibility criterion q < q v = 3 AU. 



2.3 Inclination 

Figure H shows the distribution of the cosine of the inclination for the LP comets. A spherically symmetric 
distribution would be flat in this figure, as indicated by the heavy line. Everhart (1967b) argued that 
inclination-dependent selection effects are minor. 

The inclination distribution in ecliptic coordinates is inconsistent with spherical symmetry: the \ 2 an d 
Kolmogorov-Smirnov statistics indicate probabilities of only 10 -4 and 10 -6 respectively that the distribution 
in Fig. |3|(a) distribution is flat. Some of the discrepancy may arise from the Kreutz sun-grazer group (~ 20 
members), which contributes at cosz ~ —0.8 and cos? ~ —0.4. The remaining excess at | cos i | ~ 1 may 
result from bias towards the ecliptic plane in comet searches, or from an intrinsically anisotropic distribution 
of LP comets resulting from individual stellar passages through the Oort cloud. 

The inclination distribution in Galactic c oordi nates may have a gap near zero inclination, possibly re- 



flecting the influence of the Galactic tide (§ 4.1.2 ), or confusion from the dense stellar background in the 
Galactic plane. 



2.4 Longitude of ascending node 

The distribution of longitude of the ascending node SI is plotted in Fig. ||. The flat line again indicates 
a spherically symmetric distribution. The Kreutz sun-grazers are concentrated at Q <~ 0.15 and SI ~ 4, 
and thus are responsible for the highest spike in Fig. [|. Everhart (1967a, b) concluded that O-dependent 
selection effects are likely to be negligible. The x 2 an d Kolmogorov-Smirnov statistics indicate that the 
ecliptic distribution is consistent with a flat distribution at the 80% and 90% levels respectively. 



2.5 Argument of perihelion 

Figure |^ shows the distribution of the argument of perihelion ui for the LP comets. Comets with < u> < ir 
outnumber those with ir < lu < 2ir by a factor of 395/286 = 1.38 ±0.11. This excess is partly due to the 
Kreutz group, which is concentrated at u> ~ 1.6 and Q ~ 0.15; but may also be due to observational selection 
(Everhart 1967a; Kresak 1982): comets with < lj < ir pass perihelion above the ecliptic, and are more 
easily visible to observers in the northern hemisphere. The lack of observed apparitions with u> > n reflects 
the smaller number of comet searchers in the southern hemisphere until recent times. The distribution in 
the Galactic frame has a slight excess of comets with orbits in the range sin2w > (399/282 = 0.59 ± 0.04). 
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q (AU) q (AU) 

(a) (b) 

Figure 2: Number N versus perihelion distance q for 681 LP comets, on two different scales. Data taken 
from Marsden and Williams (1993). The solid line is the estimated intrinsic distribution from Kresak and 
Pittich (1978), the dotted line is from Everhart (1967b), and the dashed line is from Shoemaker and Wolfe 
(1982). The appropriate normalizations are difficult to determine for these curves, and are chosen somewhat 
arbitrarily. 




Figure 3: The distribution of the cosine of the inclination for the 681 LP comets in (a) ecliptic coordinates, 
and (b) Galactic coordinates. A spherically symmetric distribution is indicated by the flat line. Data taken 
from Marsden and Williams (1993). 
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Figure 5: The distribution of the argument of perihelion in (a) the ecliptic frame, u, and (b) the Galactic 
frame, a), for the 681 LP comets. Data taken from Marsden and Williams (1993). 
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2.6 Aphelion direction 

Figure || shows the distribution of the aphelion directions of the LP comets in ecliptic and Galactic coordi- 
nates. 

Claims have been made for a clustering of aphelion directions around the solar antapex (e.g. Tyror 1957; 
Oja 1975), but newer analyses with improved catalogues (e.g. Lust 1984) have cast doubt on this hypothesis. 
The presence of complex selection effects, such as the uneven coverage of the sky by comet searchers, renders 
difficult the task of unambiguously determining whether or not clustering is present. Attempts to avoid 
selection effects end up subdividing the samples into subsamples of such small size as to be of dubious 
statistical value. 

Whipple (1977) has shown that it is unlikely that there are many large comet groups i.e. comets related 
through having split from the same parent body, in the observed sample though the numerous (~ 20) 
observed comet splittings makes the possibility acceptable in principle. A comet group would likely have 
spread somewhat in semimajor axis: the resulting much larger spread in orbital period P cx a 3 / 2 makes it 
unlikely that two or more members of such a split group would have passed the Sun in the 200 years for 
which good observational data exist. The Kreutz group of sun-grazing comets is the only generally accepted 
exception. 

Figures [7]a and b show histograms of comet number versus the sine of the ecliptic latitude /3 and of the 
Galactic latitude b. The ecliptic latitudes deviate only weakly from a spherically symmetric distribution and 
this deviation is likely due to the lack of southern hemisphere comet searchers. The Galactic distribution 
shows two broad peaks, centred roughly on sin b ~ ± 0.5. It will be shown that these probably reflect the 
influence of the gravitational tidal field of the Galaxy (§ 4.1.2 ), which acts most strongly when the Sun-comet 
line makes a 45° angle with the Galactic polar axis. 



2.7 Orbital elements of new comets 

For some purposes it is useful to isolate the distribution of orbital elements of the 109 new comets whose 
original semimajor axes lie in the Oort spike, x — 1/a < 10~ 4 AU _1 . The distributions of perihelion distance, 
as well as inclination, longitude of the ascending node, and argument of perihelion in Galactic coordinates 
are all shown in Fig. ^ The distribution of aphelion directions is shown in Fig. |9| 

2.8 Parametrization of the distribution of elements 

For comparison with theoretical models, we shall parametrize the observed distribution of LP comets by 
three dimensionless numbers: 

• The ratio of the number of comets in the Oort spike (1/a < 10~ 4 AU^ 1 ) to the total number of LP 
comets is denoted by "fi. This parameter measures the relative strength of the Oort spike. 

• The inverse semimajor axes of LP comets range from zero (unbound) to 0.029 AU _1 (P — 200 yr). 
Let the ratio of the number of comets in the inner half of this range (0.0145 to 0.029 AU _1 ) to the 
total be This parameter measures the prominence of the "tail" of the energy distribution. 

• Let the ratio of the number of prograde comets in the ecliptic frame to the total be 'J 3 . This parameter 
measures the isotropy of the LP comet distribution. 

We estimate these parameters using all LP comets with original orbits in Marsden and Williams (1993): 

*i = 109/289 = 0.377 ±0.042, 

* 2 = 19/289 = 0.066 ±0.016, (2) 
* 3 = 145/289 = 0.501 ±0.051. 



2 OBSERVATIONS 



9 




Figure 6: All 681 long-period comet aphelion directions on ecliptic (a) and Galactic (b) equal-area maps. 
More precisely, these are the antipodes of the perihelion directions. The crossed circle is the solar apex. 
Data taken from Marsden and Williams (1993). 
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For consistency, we based our calculation of "J 3 on the 289 comets with known original orbits, even though 
knowledge of the original orbit is not required since ^3 depends only on angular elements. If we consider all 
681 LP comets, we find VE^ = 321/681 = 0.471 ± 0.041; the two values are consistent within their error bars. 

We denote theoretical values of these parameters by \P' and compare theory and observation through the 
parameters 

Xt = ^, i= 1,2,3; (3) 
which should be unity if theory and observation agree. 



3 Theoretical background 

3.1 The Oort cloud 

The spatial distribution of comets in the Oort cloud can be deduced from the assumption that these comets 
formed in the outer planetary region and were scattered into the Oort cloud through the combined pertur- 
bations of the tide and planets (Duncan et al. 1987). These calculations suggest — in order of decreasing 
reliability — that (i) the cloud is approximately spherical; (ii) the velocity distribution of comets within the 
cloud is isotropic; in other words the phase-space distribution is uniform on the energy hypersurface, except 
perhaps at very small angular momentum where the comets are removed by planetary encounters; (iii) the 
cloud's inner edge is near 3000 AU, with a space number density of comets roughly proportional to r" 3,5 
from 3000 to 50 000 AU. 

Orbits of comets in the Oort cloud evolve mainly due to torques from the overall Galactic tidal field, 
but they are also affected by encounters with planets, passing stars and molecular clouds. Comets are also 
lost through collisions with the Sun. Through these mechanisms, between 40% (Duncan et al. 1987) and 
80% (Weissman 1985) of the original Oort cloud may have been lost over the lifetime of the Solar System, 
leaving perhaps ~ 10 12 comets (cf. eq. ^) with mass ~ 50M® (Weissman 1990) in the present-day comet 
cloud. These numbers are very uncertain. 

If the phase-space distribution of comets is uniform on the energy hypersurface, then the number of 
comets at a given semimajor axis with angular momentum less than J should be cx J 2 ; this in turn implies 
that the number of comets with perihelion in the range [q,q + dq] should be N(q)dq, where 

N(q) oc 1 - -, q < a. (4) 
a 

This distribution is modified if there are loss mechanisms that depend strongly on perihelion distance, as we 
now discuss. 

3.2 The loss cylinder 

A comet that passes through the planetary system receives a gravitational kick from the planets. The 
typical energy kick Ax depends strongly on perihelion distance (and less strongly on inclination): Aa; « 1 x 
10" 3 AU" 1 for q < 6 AU, dropping to 1 x 10" 4 AU" 1 at q ~ 10 AU and 1 x 10" 5 AU" 1 at q ~ 20 AU (Fernandez 
1981; Duncan et al. 1987). For comparison, a typical comet in the Oort spike has x < 10" 4 AU" 1 ; since 
these comets have perihelion q ~ 1 AU they receive an energy kick Aa; 3> x during perihelion passage. 
Depending on the sign of the kick, they will cither leave the planetary system on an unbound orbit, never 
to return, or be thrown onto a more tightly bound orbit whose aphelion is much smaller than the size of the 
Oort cloud. In cither case, the comet is lost from the Oort cloud. 

More generally, we can define a critical perihelion distance qg ~ 10 AU such that comets with q < qg 
suffer a typical energy kick at perihelion which is larger than the typical energy in the Oort cloud. Such 
comets are said to lie in the "loss cylinder" in phase space because they are lost from the Oort cloud within 
one orbit (the term "cylinder" is used because at a given location within the cloud, the constraint q < qi is 
is satisfied in a cylindrical region in velocity space: for highly eccentric orbits q < qg implies that the angular 
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momentum J < Jg = (2GM ( 7 ) qg) 1 ^ 2 , which in turn implies that the tangential velocity v± < Jt/r). The loss 
cylinder is refilled by torques from the Galactic tide and other sources. 

The comets in the Oort spike are inside the loss cylinder and hence must generally be on their first 
passage through the planetary system (this is why we designated the 109 comets with 1/a < 10~ 4 AU -1 
as "new" in § 2.7). The loss cylinder concept also explains why the energy spread in the Oort spike is 
much narrower than the energy spread in the Oort cloud itself: comets with smaller semimajor axes have 
a smaller moment arm and shorter period so their per-orbit angular momentum and perihelion distance 
changes are smaller; for a < 2 x 10 4 AU the perihelion cannot jump the "Jupiter barrier" i.e. cannot evolve 
from q > qi ~ 10 AU (large enough to be outside the loss cylinder) to q < 1 AU (small enough to be visible) 
in one orbital period. Thus the inner edge of the Oort spike is set by the condition that the typical change 
in angular momentum per orbit equals the size of the loss cylinder, and does not reflect the actual size of the 
Oort cloud (Hills 1981). The new comets we see come from an outer or active Oort cloud (a > 2 x 10 4 AU) 
in which the typical change in angular momentum per orbit exceeds the size of the loss cylinder. Thus, in the 
outer Oort cloud, losses from planetary perturbations do not strongly affect the phase-space distribution of 
comets near zero angular momentum (the loss cylinder is said to be "full" ) , and the equilibrium distribution 
of perihelion distances (Eq. |i|) remains approximately valid within the loss cylinder. The more massive inner 
Oort cloud (a < 2 x 10 4 AU) does not produce visible comets except during a rare comet "shower" caused 
by an unusually close stellar encounter and which perturbs them sufficiently to jump the Jupiter barrier. 
In this inner cloud, losses from planetary perturbations strongly deplete the distribution of comets at small 
perihelion distances (the loss cylinder is said to be "empty") and thus it does not contribute to the Oort 
spike. 



3.3 Energy evolution of LP comets 

Let us examine the motion of an Oort cloud comet after it enters the planetary system for the first time. 
The motion of a comet in the field of the giant planets, the Sun and the Galactic tide is quite complicated, 
but considerable analytic insight can be obtained if we make the following approximations: 

1. Depending on the sign of the energy kick from the planets, the comet will either be ejected from the 
Solar System or perturbed onto a much more tightly bound orbit (a ~ 10 3 AU). In either case the 
Galactic tide plays no further significant role, and can be neglected. 

2. The influence of the planets is concentrated near perihelion, where the moment arm of the comet 
is small. Thus the angular momentum, orbit orientation, and perihelion distance are approximately 
conserved during perihelion passage; the only significant change is in the orbital energy. More precisely, 
the typical changes in angular elements caused by a planet of mass M p and semimajor axis a p to the 
orbit of a comet with perihelion q < a p are 

Ai ~ AS1 ~ Auj ~ Aq/a p ~ M p /M Q , (5) 

while the fractional change in energy is Ax/x ~ [aM p / cl p Mq), which is much larger for an Oort cloud 
comet (by a factor a/a p ~ 10 3 -10 4 ). 

3. The orbit of any LP comet looks nearly parabolic when it passes through the planetary system, so the 
distribution of energy changes during perihelion passage is approximately independent of the comet's 
orbital energy. Therefore we may define a function p( Ax)dAx, the probability that the energy change 
per perihelion passage due to planetary perturbations is in the interval [Ax, Ax + dAx]. The function 
p(Ax) is an implicit function of the inclination, perihelion distance, argument of perihelion, etc. but 
as we have seen these change much more slowly than x and so can be considered constants. The 
properties of p(Ax) are discussed by Everhart (1968); p is approximately an even function of Ax [the 
odd component is smaller by O(M p /M )] and as |Ax| — > oo, p(Ax) oc |Acc| -3 — although despite this 
extended tail p is often approximated by a Gaussian. If q < a p then the typical energy change due to a 
single planet is (|Ax|) ~ (M p /Mo)a~ 1 . Plots of the second moment oip(Ax), averaged over argument 
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of perihelion, as a function of perihelion distance and inclination are given by Fernandez (1981) and 
Duncan et al. (1987). 

4. Comets on escape orbits (x < 0) are lost from the Solar System. The appropriate boundary condition 
at large x is less clear. Our other approximations fail when x becomes comparable to the inverse 
semimajor axes of the planets, which occurs when they become SP comets at x = x sp — 0.029 AU -1 . 
SP comets will continue to random walk in energy — some becoming LP comets once again — but the 
other orbital elements will also evolve at a comparable rate, so the approximation of a one-dimensional 
random walk is no longer valid. Fortunately we shall find that the fraction of new comets that survive 
fading and ejection to become SP comets is small enough that the details of their evolution are unlikely 
to affect the overall distribution of LP comets (cf. Tables [j] and ; for most purposes we can simply 
assume that LP comets reaching x sp are lost. 

5. The orbital periods of most LP comets are sufficiently long that the orbital phases of the planets and 
hence the energy kicks that the comet receives from them are uncorrelated at successive encounters. 
Thus the evolution of the comet energy can be regarded as a Markov process or random walk.Q 

6. A comet may fade as its inventory of volatiles is depleted or may be disrupted by various mechanisms. 
We shall use the term "fading" to denote any change in the intrinsic properties of the comet that would 
cause it to disappear from the observed sample. We parametrize this process by a function <3>. m £ [0, 1], 
m = 1, 2, . . . ($i = 1), the probability that a visible new comet survives fading for at least m perihelion 
passages. There are two closely related functions: the probability that the comet survives fading for 
precisely m perihelion passages, 

<t> m = $ m - $ m+ i; (9) 

and the conditional probability that a comet that survives m passages will fade before the (m -I- l) st 
passage, 

tPm. = "f— = 1 ^ ■ (10) 

With these approximations the evolution of the energy of LP comets can be treated as a one-dimensional 
random walk. We assume that visible new comets arrive directly from the Oort cloud with original energy 
x = 0, at a rate N per year. Let f m {x)dx be the number of visible LP comets per year with original energy 
in the range [x,x + dx] which are returning on their m th perihelion passage (thus fi(x) = NS(x), where 8(x) 
is a delta function). Then in a steady state we must have 

/>oo 

fm+ifa) = (1 - tym) I fm(y)p(x - y)dy, m > 1, x > 0, (11) 
Jo 

3 This argument can be made more precise (Chirikov and Vecheslavov 1989; Pctrosky 1986; Sagdeev and Zaslavsky 1987). 
Assume that the energy kick received by a comet is proportional to the sine of the orbital phase of Jupiter relative to the comet. 
Let x n be the original energy (in AU — 1 ) just before the n perihelion passage and let g n be the orbital phase of Jupiter at 
that passage. Then 

2-7T —3/2 

x n+1 = x n + ei 3ing„, g n +i = 9n + — x n+ \ , (6) 

where Pj is Jupiter's orbital period in years and e = 2 — 1 ' 8 ei is the rms energy change for a comet with small perihelion. 
Writing x n = xo + Sx n where 8x n is small, the map becomes 

In+i = In + £i/sing„, g n +i = g n — I n +i, (7) 

where I n = fSx n and / = 3n / ' (Pjx^J 2 ). This is the standard map, which exhibits global chaos when |ei/| -Sj 1; this in turn 
implies that the energy kicks received by the comet are effectively random. The condition |ei/| ^ 1 can be re-written as 

a > 2-^4/5(3^-2/5 = 20 AU ( 5 X W ~* AU "' \ ^ , (8) 

where e RJ 5 X 10 — 4 AU — 1 is a typical energy kick for comets with perihelia inside Jupiter's orbit (Fernandez 1981; Duncan 
et al. 1987). 

Thus the orbits of all LP comets with perihelion q < aj are expected to be chaotic. 
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which can be solved successively for f2(x), fd(x), ■ ■ ■• The total number of LP comets with energies in the 
interval [x, x + dx] is f(x)dx, where f(x) = X)m=o theoretical predictions of f{x) are to be compared 

with the observed distribution of LP comets in Fig. [j]. 

The simplest version of this problem is obtained by assuming that there is no fading (ip m = 0) and 
that the energy changes by discrete steps ±e with equal probability [p(x) — + e) + \5{x — e)]. In this 
case the possible values of the energy are restricted to a lattice x = (j — l)e, where j is an integer, and 
the random walk is identical to the gambler's ruin problem (Kannan 1979; Feller 1968). The end-state of 
ejection (j = 0) corresponds to bankruptcy; if in addition we assume that there is an absorbing boundary at 
x sp = (j sp — l)e, then evolving to an SP comet corresponds to breaking the house. Thus, for example, the 
probabilities that an LP comet with energy (j — l)e will be eventually be ejected or become a short-period 
comet are respectively 

Pej = 1 - — , Psp = — , (12) 

Jsp Jsp 

and the mean number of orbits that the comet will survive is 

(m) =j(j sp -j), (13) 

A new comet has j — 1 and its mean lifetime is therefore (m) = j sp — 1; the ratio of new to all LP comets 
observed in a fixed time interval is 

** = -!- = — - — . (14) 
1 (m) j sp -l 1 ' 

There are also explicit expressions for the probability that the comet is ejected or becomes an SP comet at 
the m th perihelion passage (Feller 1968). 

The gambler's ruin problem is particularly simple if there is no boundary condition at large x (x sp — ► oo), 



which is reasonable since few comets reach short-period orbits anyway (§ 5.2.1). The probability that a new 
comet will survive for precisely m orbits is then 



Poj(m) = 2^Um+iJ' m ° dd ' 

=0 to even; (15) 

for to 3> 1, p c j(to) — > (2/7r) 1 / 2 m -3 / 2 for to odd, and zero otherwise. The mean lifetime X)m=i m Pcj(™) is 
infinite, and the probability that a comet will survive for at least m orbits is ~ mT 1 / 2 for large m. 

When using the gambler's ruin to model the evolution of LP comets, we take e ~ 5 x 1(T 4 All" 1 , which 
is the rms energy change for comets with perihelion between 5 and 10 AU (Fernandez 1981; Duncan et al. 



1987), and z sp = 0.029 AIT 1 (P = 200 yr); thus j sp ~ 60. Eq.yjthen predicts V\ = 0.017; the ratio of the 
predicted to the observed value for this parameter (cf. Eq. ||) is 

X t = -± = 0.051 ± 0.006. (16) 

The gambler's ruin model predicts far too few comets in the Oort spike relative to the total number of LP 
comets. 

This simple model also makes useful predictions about the inclination distribution of LP comets. The 
distribution of new comets is approximately isotropic, so there are equal numbers of prograde and retrograde 
new comets. Since prograde comets have longer encounter times with the planets, they tend to have larger 
energy changes than retrograde comets. Equation |lj predicts that the ratio of prograde to retrograde LP 
comets should be roughly the ratio of the rms energy change for these two types, e re tro/epro —2—3. The 
fraction of prograde comets should then be ^3 = 1/(1 + e pro /e re tro) — 0.3. The ratio of the predicted to the 
observed value for this parameter (cf. Eq. ||) is 

A 3 = = 0.58 ± 0.06. (17) 



3 THEORETICAL BACKGROUND 



14 



The gambler's ruin model predicts too few prograde comets. 

More accurate investigations of this one-dimensional random walk have been carried out by many authors. 
Oort (1950) approximated p(Ax) by a Gaussian and assumed tp m = k = constant and found a good fit to 
most of the energy distribution for k — 0.014; however, he found that the number of new comets was larger 
than the model predicted by a factor of five, and hence was forced to assume that only one in five new comets 
survive to the second perihelion passage — in other words ipi = 0.8, tpm = 0.014 for to > 1. Kendall (1961) 
and Yabushita (1979) have analyzed the case x sp —* oo, ip m = k = constant, p(y) oc exp (— 2 1 / 2 |?/|/er) , where 
a is the rms energy change per perihelion passage. In this case Eq. [ll] can be solved analytically to yieldQ 

Ol/ 2 AT 

f{x) = N5(x) + (1 - fc 1/2 ) exp [ - (2k) 1 / 2 x/a] ; (18) 

using this result Kendall derives a reasonable fit to the data if k = 0.04 and one in four to six new comets 
survive to the second perihelion — results roughly compatible with Oort's. This model predicts a ratio of new 
comets to all LP comets observed in a fixed time interval given by 

*i = t#vT = kl/2 - ^ 
J f{x)dx 

Yabushita (1979) gave analytic formulae for p e j(m) for this model, and showed that the probability that a 
comet will survive for at least to orbits is ~ exp(— krnj/m 1 ^ 2 for large m. Whipple (1962) examined survival 
laws of the form <p m oc m~ a (the proportionality constant is determined by the condition that ^ m <p m = 1) 
and found a good fit to the observed energy distribution with a ~ 1.7. Everhart (1979) used a distribution 
p(y) derived from his numerical experiments and found <j> m ~ 0.2 for all to > 1; in other words only one in 
five comets survived to the second perihelion passage but the fading after that time was negligible. 

For some purposes the random walk can be approximated as a diffusion process; in this case the rele- 
vant equations and their solutions are discussed by Yabushita (1980). Bailey (1984) examines solutions of 
a diffusion equation in two dimensions (energy and angular momentum) and includes a fading probability 
that depends on energy rather than perihelion number — which is less well-motivated but makes the equa- 
tions easier to solve (he justifies his fading function with an a posteriori "thermal shock" model, in which 
comets with large aphelia are more susceptible to disruption because they approach perihelion with a lower 
temperature) . Bailey finds a good fit to the observed energy distribution if the fading probability per orbit 
is 

cj){x) = 0.3[1 + (z/0.004 AU) 2 ]~ 3/2 . (20) 

Emel'yanenko and Bailey (1996) have modeled the distribution of LP comets using a Monte Carlo model 
with ip m = k = constant plus an additional probability per orbit k* that the comet is rejuvenated. Their 
preferred values are k = 0.3 and k* = 0.0005. 

The most complete model of LP comet evolution based on a random walk in energy is due to Weiss- 
man (1978, 1979, 1980). His Monte Carlo model included the gravitational influence of the planets, non- 
gravitational forces, forces from passing stars, tidal disruption by the Sun, fading and splitting. In his 
preferred model, 15% of the comets have zero fading probability, and the rest had a fading probability of 0.1 
per orbit. At this cost of this somewhat ad hoc assumption, Weissman was able to successfully reproduce 
the semimajor axis, inclination, and perihelion distributions. 

The one-dimensional random walk is a valuable tool for understanding the distribution of LP comets. 
However, some of its assumptions are not well-justified: (i) Secular changes in perihelion distance, argument 
of perihelion, and inclination at each perihelion passage accumulate over many orbits and can lead to sub- 
stantial evolution of the orientation and perihelion (Quinn et al. 1990; Bailey et al. 1992; Thomas and 
Morbidelli 1996). (ii) Although the probability distribution of energy changes p(y) is approximately an even 
function [(y 2 ) 1 ^ 2 is larger than (y) by 0(M p /M@)], the random changes in energy due to the second moment 
grow only as to 1 / 2 where to is the number of orbits, while the systematic changes due to the first moment 
grow as to. Thus the small asymmetry in p(y) may have important consequences. 

4 In the limit of zero disruption, k = 0, Eq. ^if yields f(x) = N[S{x) + 2 1 / 2 /a] ; in other words, the energy distribution of 
observed LP comets should be flat, apart from the Oort spike. 
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3.4 The fading problem 

All the investigations described in the previous subsection reach the same conclusion: if the LP comets are 
in a steady state then we cannot match the observed energy distribution without unexpectedly strong fading 
after the first perihelion passage. Therefore either (i) the comet distribution is not in a steady state, which 
almost certainly requires rejecting most of the Oort model^], or (ii) we must postulate ad hoc fading laws 
and abandon the use of the energy distribution as a convincing test of the Oort model. This is the fading 
problem. 

Fading can arise from many possible mechanisms but the most natural hypothesis is that the comet's 
brightness fades sharply because its inventory of volatiles is depleted during the first perihelion passage. 
Oort and Schmidt (1951) have argued that this hypothesis is supported by the observation that new comets 
have strong continuum spectra due to dust entrained by the gases from a volatile component, and that the 
decline of brightness with increasing heliocentric distance is much slower for new comets. Many authors have 
looked for evidence that new comets differ in composition or brightness from older LP comets, with mixed 
results; Whipple (1991) summarizes these investigations by saying that the Oort-Schmidt effect is "fairly 
well confirmed" . 

Fading is much slower after the first perihelion passage, as exemplified by the long history of Hallcy's 
comet. Whipple (1992) concludes that there is no strong evidence that older (i.e. shorter period) LP comets 
have faded relative to younger LP comets, consistent with theoretical estimates that 10 3 -10 4 orbits are 
required for moderate-sized comets to lose their volatiles (Weissman 1980) and the lack of strong systematic 
trends in the brightness of SP comets. 

Comets may also fade if they disrupt or split. After splitting, the fragments are fainter and hence less 
likely to be visible, and in addition lose their volatiles more rapidly. Moreover, young comets are more likely 
to split than old ones: Weissman (1980) gives splitting probabilities per perihelion passage of 0.10 ± 0.04 for 
new comets but only 0.045 ± 0.011 for LP comets in general. The cause of splitting is not well understood, 
except in some cases where splitting is due to tidal forces from a close encounter with a giant planet. 

Finally, we note that LP comets are responsible for 10-30% of the crater production by impact on Earth 
(Shoemaker 1983). The observed cratering rate can therefore — in principle — constrain the total population 
of LP comets, whether or not they have faded; however, this constraint is difficult to evaluate, in part because 
estimates of comet masses are quite uncertain. 

4 Algorithm 

We represent each comet by a massless test particle and neglect interactions between comets. The orbit of 
the test particle is followed in the combined gravitational fields of the Sun, the four giant planets, and the 
Galactic tide. We assume that the planets travel around the Sun in circular, coplanar orbits. We neglect the 
terrestrial planets, Pluto, the small free inclinations and eccentricities of the giant planets, and their mutual 
perturbations as there is no reason to expect that these play significant roles in the evolution of LP comets. 

4.1 Equations of motion 

The equation of motion of the comet can be written as 

r = F© + Fplancts + Ftido + F thcr, (21) 

where the terms on the right side represent the force per unit mass from the Sun, the planets, the Galactic 
tide, and other sources (e.g. non-gravitational forces). 

5 There are advocates of this position (see Bailey 1984 for references), but we are not among them. 
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4.1.1 The planets 

We shall employ two frames of reference: the barycentric frame, whose origin is the center of mass of the 
Sun and the four planets, and the heliocentric frame, whose origin is the Sun. In the barycentric frame, 

GM n , , ^ GM, 



Wts = - r «) - E rr-^ r ~ r ^' ( 22 ) 

p 



where r, r Q , and r p are the positions of the comet, the Sun, and planet p. In the heliocentric frame, the Sun 
is at the origin and 

„ „ GM Q x - GM V , . x - GMr, 

F Q + F plancts = —j-j^r - J2 jT"^ ( r - r P ) ~ E ( 23 ) 

the last sum is the "indirect term" that arises because the heliocentric frame is not inertial. 

The heliocentric frame is useful for integrating orbits at small radii, |r| < |r p |, because it ensures that the 



primary force center, the Sun, is precisely at the origin (see §4.1.4). It is not well-suited for integrating orbits 
at large radii, |r| > r p |, because the indirect term does not approach zero at large radii, and oscillates with 
a period equal to the planetary orbital period — thereby forcing the integrator to use a very small timestep. 
In the integrations we switch from heliocentric to barycentric coordinates when the comet radius |r| exceeds 
a transition radius r c ; tests show that the integrations are most efficient when r c = 10 AU. 

The code tracks close encounters and collisions between comets and planets. A close encounter with a 
planet is defined to be a passage through a planet's sphere of influence 



M, 



2/5 



where a p is the planet's semimajor axis. Each inward crossing of the sphere of influence is counted as 
one encounter, even if there are multiple pericenter passages while the comet remains within the sphere of 
influence. A close encounter with the Sun is defined to be a passage within 10 solar radii. 



4.1.2 The Galactic tide 

The effects of the Galactic tide on comet orbits are discussed by Heisler and Tremaine (1986), Morris and 
Mullcr (1986), Torbett (1986), and Matese and Whitman (1989). Consider a rotating set of orthonormal 
vectors {ej,e^,ej}. Let ej point away from the Galactic center, e^ in the direction of Galactic rotation, 
and e^ towards the South Galactic Pole (South is chosen so that the coordinate system is right-handed). 
The force per unit mass from the tide is (Heisler and Tremaine 1986) 

F tido = {A- B)(3A + B)xe s -(A- B) 2 ye.y - [inGp - 2{B 2 - A 2 )]~ze~ z , (25) 

where po is the mass density in the solar neighborhood, and A and B are the Oort constants. We take 
A = 14.4 ± 1.2 km s" 1 kpc" 1 and B = -12.0 ± 2.8 km s" 1 kpc" 1 (Kerr and Lynden-Bell 1986). The 
local mass density is less well-known. Visible matter (stars and gas) contributes about 0.1 M pc~ 3 , but 
the amount of dark matter present in the solar neighbourhood remains controversial. If the dark matter is 
distributed like the visible matter, then the dark/visible mass ratio P is between and 2 (Oort 1960; Bahcall 
1984; Kuijken and Gilmore 1989; Kuijken 1991; Bahcall et al. 1992). We adopt p a = 0.15 M Q pc~ 3 in 
this paper, corresponding to P = 0.5. 

With these values of A, B and po, the AirGpo term of Eq. ^ exceeds the others by more than a factor of 
ten, and from now on we shall neglect these other terms. The dominant component of the tidal force arises 
from a gravitational potential of the form 



Vtidc = 2ttGp z 2 . 



(26) 
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In practice, of course, the local density po varies as the Sun travels up and down, in and out, and through 
spiral arms during its orbit around the Galaxy. The amplitude of this variation depends strongly on the 
unknown distribution and total amount of disk dark matter. The maximum-to-minimum density variation 
could be as large as 3:1 (Matese et al. 1995) but is probably considerably smaller, with a period around 
30 Myr [close to ^(n/Gpo) 1 ^ 2 , the half-period for oscillations in the potential ( p6[ ) ]. We are justified in 
neglecting these variations in poi because the typical lifetime of LP comets after their first apparition is only 
1.4 Myr (see Table III below), which is much shorter. 



4.1.3 Encounters with stars and molecular clouds 

Our model neglects the effects of passing stars on LP comets, for two main reasons: (i) The delivery rate 
of Oort cloud comets to the planetary system due to Galactic tides is higher than the rate due to stellar 
encounters by a factor 1.5-2 (Heisler and Tremaine 1986; Torbett 1986), except during rare comet showers 
caused by an unusually close passage, during which the delivery rate may be enhanced by a factor of twenty 
or so (Hills 1981; Heisler 1990); we feel justified in neglecting the possibility of a comet shower because they 
only last about 2% of the time (Heisler 1990); (ii) The effects of stellar encounters are highly time- variable 
whereas the strength of the tide is approximately constant over the typical lifetime of LP comets; thus by 
concentrating on the effects of the tide we focus on a deterministic problem, whose results are easier to 
interpret. 

The effects of rare encounters with molecular clouds are highly time- variable, and difficult to estimate 
reliably because the properties of molecular clouds are poorly known (Bailey 1983; Drapatz and Zinnecker 
1984; Hut and Tremaine 1985; Torbett 1986). Therefore we shall also assume that the present distribution 
of LP comets has not been affected by a recent encounter with a molecular cloud. 



4.1.4 Regularisation 

Integrating the orbits of LP comets is a challenging numerical problem, because of the wide range of timescales 
(the orbital period can be several Myr but perihelion passage occurs over a timescale as short as months) and 
because it is important to avoid any secular drift in energy or angular momentum due to numerical errors. 
We have used the Kustaanheimo-Stiefel (K-S) transformation to convert Cartesian coordinates to regularized 
coordinates and have carried out all of our integrations in the regularized coordinates. A requirement of 
K-S regularisation is that the frame origin must coincide with the primary force centre, which is why we use 
heliocentric coordinates at small radii. 

The numerical integrations were carried out using the Bulirsch-Stoer method, which was checked using 
a fourth-order Runge-Kutta-Fehlberg algorithm. All integrations were done in double-precision arithmetic. 



4.2 Non-gravitational forces 

The asymmetric sublimation of cometary volatiles results in a net acceleration of the nucleus. These non- 
gravitational^] (NG) forces are limited to times of significant outgassing (i.e. coma production), and remain 
small even then. 

Non-gravitational forces are difficult to model. Their strength obviously depends on the comet's distance 
from the Sun, but displays less regular variability as well: gas production may vary by a factor of 2 or more 
between the pre- and post-perihelion legs of the orbit (Sekanina 1964; Festou 1986), and jets and streamers 
are observed to evolve on time scales of less than a day (Festou et al. 1993b), suggesting that NG forces 
change on similar time scales. Further complications arise from the rotation of the nucleus, which is difficult 
to measure through the coma, and which may be complicated by precession (Wilhelm 1987). 

6 Traditionally, the term "non-gravitational forces" has been reserved for the reaction forces resulting from the uneven 
sublimation of cometary volatiles, and it will be used here in that manner. Other factors of a non-gravitational nature have been 
considered, including radiation and solar wind pressure, drag from the interplanetary /interstellar medium, and the heliopausc, 
but were found to be negligible in comparison to the outgassing forces (Wiegert 1996). 
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The NG acceleration Fj ct is written as 

Fjot = ^ie 1 + F 2 e 2 + F 3 e 3 , (27) 

where ei points radially outward from the Sun, e 2 lies in the orbital plane, pointing in the direction of orbital 
motion and normal to ei, and e 3 = ei x e 2 . A naive model of NG accelerations, which is all the data allows, 
assumes that the short timescale components are uncorrelated and cancel out, leaving only fairly regular, 
longer timescale components as dynamically important. We shall use the Style II model of Marsden et al. 
(1973), which assumes that accelerations are symmetric about perihelion, and can be represented by 

Fxir) = A l9 {r), F 2 (r) = A 2 g(r), F 3 (r) = A 3 g(r). (28) 

Here {Ai, A 2 , A3} are independent constants, and g(r) is a non-negative function describing the dependence 
on the comet-Sun distance r. The form of g(r) is based on an empirical water sublimation curve by Delscmmc 
and Miller (1971), 



'"0 



1+1 - 



(29) 



where m = 2.15, k = 4.6142, n = 5.093, r Q = 2.808 AU and a is chosen to be 0.1113 so that g(l AU) = 1. 
Note that g(r) is roughly proportional to r~ m w r~ 2 for r C ro. At r 3> r$, g(r) drops much faster than 
the simple inverse square that describes the incident solar flux. 

The constants Ai are determined by fitting individual comet orbits (Marsden et al. 1973); the value of 
Ai is typically 10" 7 to 10~ 9 AU day -2 , \A 2 \ is typically only 10% of \Ai\, and A3 is consistent with zero. 

4.3 Initial conditions 

4.3.1 Initial phase-space distribution 

The distribution of comets in the Oort cloud is only poorly known, although it is plausible to assume that 
the cloud is roughly spherical and that the comets are uniformly d istri buted on the energy hypersurface in 




phase space, except possibly at very small angular momenta (cf. §3T). Then the phase-space density is a 
function only of L = (GAf a) 1 / 2 , which we assume to be 

L<L- = (GMgfl-) 1 / 2 , 

L_ < L < L+, (30) 
L> L+ = (GM Q a + ) 1 / 2 , 

where fo and a are constants, and a_ and a+ are the inner and outer edges of the Oort cloud, respectively. 
We show below (footnote 7) that the total number of Oort cloud comets with semimajor axes in the range 
specified by [L, L + dL] is (2tt) 3 f(L)L 2 dL; this in turn implies that the number density of comets is oc r a 
for a_ «C r <C a+. 

Simulations of the formation of the Oort cloud by Duncan et al. (1987) suggest that the number density 
of Oort cloud comets is cx r -( 3 - 5±a5 ) between 3000 and 50 000 AU. Thus we set a = -3.5, a- = 10 000 AU 
and a+ = 50 000 AU. The inner edge of the cloud was placed at 10 000 AU instead of 3000 AU because 
comets with a < 10 000 AU cannot become visible except in occasional comet showers, yet would consume 
most of the computer time in our simulation. 

If the comets are uniformly distributed on the energy hypersurface, the fraction of cloud comets with 
perihelion less than q < a is J 2 (q)/L 2 = 2q/a = 0. 003(^/40 AU)(25, 000 AU/a) (which is consistent with 
J N(q)dq as given by Eq. |j). Since the effects of the planets decline rapidly to zero when q > 40 AU, 
only a small fraction of cloud comets are influenced by planetary perturbations. Therefore to avoid wasting 
computer time we analyze the motion of comets with larger perihelion distance analytically, as we now 
describe. 
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4.3.2 Orbit-averaged evolution 

For comets in the Oort cloud, the tidal potential (|26| ) is much smaller than the Kepler Hamiltonian Hkcp = 
— ^GMqI a. Thus the evolution of the comet under the Hamiltonian i?Kop + ^tidc can be approximately 
described by averaging Vtide over one period of a Kepler orbit to obtain the orbit-averaged Hamiltonian 
(Hcisler and Tremaine 1986) 

GM 

H m = — ^ + vrGpo a 2 sin 2 ? (1 - e 2 + 5e 2 sin 2 u>); (31) 

2a 

here i and u) are the inclination and argument of perihelion measured in the Galactic frame. It is useful to 
introduce canonical momenta 

L = (GM a) 1/2 , J = [GM Q a(l - e 2 )] 1/2 , J s = J cos I (32) 

and their conjugate coordinates 

/, u>, a (33) 

Here J is the usual angular momentum per unit mass, Jj is its component normal to the Galactic plane, / 
is the true anomaly and tl is the longitude of the ascending node on the Galactic planeQ In terms of the 
canonical coordinates and momenta the orbit-averaged Hamiltonian is 

(GMq) 2 irp L 2 2 2 r 2 2 7 2\ • 2-1 ,n A \ 

2lX~ GMjl 2 ^ ~ J i>V J + 5 ( L -JjsmwJ. (34) 

The canonical variables / and f2 are absent from Eq. [m], so the conjugate momenta L and Ji are conserved. 
The conservation of L implies that semimajor axis is conserved as well. The solution of the equations of 
motion ( |3^ ) is discussed by Heisler and Tremaine (1986) and Matese and Whitman (1989) but is not needed 
for our purposes. 

The rate of change of angular momentum is given by 

J = =77-, (35) 



571-po L 2 



GM 2 J 2 



{J 2 - Ji){L z - J ) sin2J), (36) 



~ GM% 



e 2 J L 4 sin 2 ?sin2w, (37) 



We now define the "entrance surface" to be the boundary of the region of phase space with J < Je(o-) = 
[2GM Q g£(ft)] 1 / 2 . We shall follow cometary orbits only after they cross the entrance surface. We choose 
cje = max(gi,g2) where qi and (72 reflect two criteria that must be satisfied by the entrance surface: (1) 
Planetary perturbations must be negligible outside the entrance surface; we take q± — 60 AU since outside 
this perihelion distance the rms fractional energy change per orbit caused by the planets is < 0.1% for a 
typical Oort cloud comet. (2) The orbit-averaged approximation for the effects of the Galactic tide must 
be reasonably accurate outside the entrance surface; thus we demand that Je must exceed 77 > 1 times the 
maximum change in angular momentum per orbit, which in turn requires 

ML)>V^L\ or a 2 = V 2 ^a\ (38) 

where we have assumed e ~ 1. In this paper we take 77 = 3. 

7 At this point we may prove a result mentioned in § [4.3. l| : if the phase-space density is / = f(L) then the total number of 
comets in the range [L, L + dL] is dN = f(L)dL dj dj s J 2 " dQ J 2 " dw J* 2 " df = (2tt) 3 f(L)L 2 dL. 
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The semimajor axis a\ t 2 where q\ = is 

- 2 , lxl 0.A U (|)" 7 (^_) ,/7 (^_ J )- 1 ' 7 . (40, 

Thus 

60 AU where a < ai^ 

8-= 4 C « ATT^ « V U . (41) 



60 AU where a > a\ 2 

>1,2/ 



4.3.3 The flux of comets into the entrance surface 



We have assumed in §4.3.1 that the phase-space density is a function only of energy or semimajor axis, 
/ = f(L). This assumption is not in general correct for small angular momentum, where the comets 
are removed by planetary encounters. However, all we require is the flux into the entrance surface, most of 
which arises from comets whose angular momentum is steadily decreasing under the influence of the Galactic 
tide. Such comets are unaffected by the planets until after they cross the entrance surface, and hence the 
assumption that / = f(L) should be approximately correct. 

Let Je, Jz,&,&, /) dL dJz d(l doj df be the current of Oort cloud comets crossing into the entrance 
surface Je at a given point. Then from Eq. |3^ 

A/r T t c\ ~ £\ I —Jf(L), where J < 

f) = { i otherwisc 

-§j^^f(L)(J 2 E-Ji)(L 2 -J 2 E )sm2Cu where sin 2Co > 0, 

otherwise. 

(43) 

In our simulations, the initial orbital elements of the comets are drawn from the distribution described 
by <£>, using the energy distribution (Eq. f30|). 



4.4 End-states 

End-states may represent the loss or destruction of a comet or simply an intermediate stopping point, from 
which the simulation can subsequently be restarted. The possible end-states are: 

Collision The distance between the comet and the Sun or one of the giant planets is less than that object's 
physical radius. To ensure that we detect collisions, when a comet is close to a Solar System body we 
interpolate between timesteps using a Keplerian orbit around that body. 

Ejection The comet is either (i) leaving the Solar System on an orbit which is unbound i.e. parabolic or 
hyperbolic with respect to the Solar System's barycentre, or (ii) has ventured beyond the last closed 
Hill surface around the Sun, and is thus considered stripped from the Solar System by the action of 
passing stars, molecular clouds, etc. In either case, the simulation is not terminated until the comet 
is at least 10 5 AU from the Sun, to allow for the possibility that subsequent perturbations will result 
in the comet losing energy and returning to a "bound" state. 

Exceeded age of Solar System The elapsed time has exceeded the age of the Solar System, 5 x 10 9 yr. 
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Exceeded orbit limit The comet has completed more than 5000 orbits without reaching one of the other 
end-states. The integration is terminated and the orbital elements are saved for later examination. 
This is a safeguard to prevent extremely long-lived comets from consuming excessive computer time. 

Faded The comet is considered to have faded through loss of volatiles, splitting or other mechanisms, and 
is no longer bright enough to be observed, even if its orbit should carry it close to the Sun or Earth. 
We shall investigate various empirical models for fading. The fading end-state is not activated in any 
simulations unless explicitly mentioned in the accompanying text. 

Perihelion too large The comet's perihelion q has evolved beyond some limit, usually taken to be 40 AU, 
and is moving outwards under the influence of the tide. Such a comet is unlikely to become visible in 
the near future. 

Short-Period The comet's orbital period has decreased below 200 yr: it has become a short-period comet. 
Continued planetary perturbations may cause short-period comets to evolve back into LP comets, but 
we shall see that the fraction of comets that reach this end-state is very small (at most a few percent; 
see Tables [| and [II), so former short-period comets are not a significant contaminant. 



Visible The comet passes within 3 AU of the Sun for the first time, an event we shall call the first apparition. 
Such comets continue to evolve, but the first apparition provides a useful intermediate stopping point 
for the simulations. 



5 Results 

We follow the trajectories of our sample comets from the time they cross the entrance surface until they 



reach one of the end-states in §4.4. We divide the evolution into two stages: the pre- visibility stage, which 



lasts until the comet first becomes visible, that is, until their first passage within 3 AU of the Sun (the first 



apparition; cf. §2.2); and the post-visibility stage, which lasts from the first apparition until the comet 
reaches one of the other end-states. 

We call the set of LP comets at their first apparition the V\ comets. Similarly, those making their 771 th 
apparition are called the V m comets. The union of the sets of orbital elements Vi, V2, ■ ■ ■ is called the 
comets. 

We intend to compare the distribution of elements of the Vy comets to the observed distribution of 
elements of new comets, and the Voo comets to the visible LP comets. Note that the Voo comets represent 
all apparitions of a set of Oort cloud comets that first crossed the entrance cylinder in a given time interval, 
while the observations yield all the comets passing perihelion in a given time interval — one is a fixed interval 
of origin and the other is a fixed interval of observation. However, in a steady state these two distributions 
must be the same except for normalization. 

For some purposes it is useful to estimate this normalization, i.e., to estimate the time interval to which 
our simulation corresponds. To do this, we first estimate the number of perihelion passages per year of new 
comets with q < q v = 3 AU, which we call $ ncw . Kresak and Pittich (1978) find that the rate of long-period 
comets passing within Jupiter's orbit (5.2 AU) is 25 yr • Taking a round number of 10 yr" 1 passing within 
3 AU and assuming one in three of these is new (Festou et al. 1993b), we find $ ncw ~ 3 yr -1 . The number 
of V\ comets produced in our simulation (see below) is 1368; hence our simulation corresponds to a time 
interval 

'3 yr" 1 " 



t. = 450 yr . (44) 

\ ^ncw / 

The total number of comets crossing the entrance surface in our simulation is 125 495. Using our assumed 
form for the semimajor axis distribution of comets in the Oort cloud (Eq. p0|, with a = —3.5), and our 



formula for the flux through the entrance cylinder (Eq. 43 ) , we may deduce that the normalization constant 
in Eq. ^o| is fo = 2.3 x 10 12 (<I> n ew/3 yr -1 ) and the total population of the Oort cloud is 

iVoort - 5.1 x 10 n ($ ncw /3 yr" 1 ), (45) 
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from 10,000 AU to 50,000 AU. Extrapolating in to 3000 AU yields a population 5% higher. For comparison, 
Heisler (1990) found that 0.2 new comets per year with perihelion < 2 AU are expected per 10 11 Oort cloud 
objects outside 3000 AU; this corresponds to an Oort cloud population that is a factor of two higher than 
the one given in Eq. ^5|. Of course, these estimates depend strongly on uncertain assumptions about the 
extent of the inner Oort cloud. 



5.1 Pre- visibility evolution 

The dynamically new or Vi comets can be used as a starting point for any investigation of phenomena that 
only affect the comet after its first apparition (non-gravitational forces, fading, etc.). The elements of the 
V\ comets are measured in the barycentric frame 200 AU from the Sun. 

The simulations reported here followed the evolution of 125 495 Oort cloud comets that crossed the 



entrance surface. The orbital elements at the entrance surface were determined as described in §4.3. Of 
the comets crossing the entrance surface, 84% had minimum perihelion distances (determined from contours 
of the averaged Hamiltonian in Eq. |34|) greater than 40 AU, too far outside the planetary system to suffer 
significant (> 1%) perturbations in scmimajor axis from the planets. These comets were transferred to the 
Perihelion too large end-state. The orbits of the remaining 20 286 comets were followed in the field of 
the Galactic tide and the Sun and planets. Table | shows the distribution of these comets among the various 
end-states; 1368 or 6.7% became V\ comets. Only 57 comets triggered the Exceeded Orbit Limit flag 



(see § L4), set at 5000 revolutions; these are discussed further in § |5.1.l| . These computations consumed 



eight weeks of CPU time on a 200 MHz workstation. 

During the pre-visibility stage there were 729 close encounters (Eq. |24| ) with the giant planets by 343 
individual comets, distributed as shown in Table 

A scatter plot of perihelion distance versus original semimajor axis for the V\ comets is shown in Fig. |l(i|a. 
There is a sharp lower bound to the distribution of semimajor axes for comets with perihelion q < 1 AU, 
which is due to the Jupiter barrier ( |3.2| ). This lower bound shifts to smaller semimajor axes at larger 
perihelion distances, since the angular momentum "hop" over the Jupiter barrier is smaller. As a result the 
number of V\ comets as a function of perihelion distance (Fig. |i~0|b) is approximately flat, as predicted by 
Eq. [|, but slightly larger for large perihelion distance. In comparison, the distribution of perihelion distances 
for the observed new comets (Fig. ||a) is not flat, but this is probably a result of the strong selection effects 
acting against comets with large perihelion. 

The distribution of original semimajor axes of the V\ comets is shown in Fig. |TT|a. The cutoff at l/a = 
2 x 10~ 5 AU -1 or a = 50 000 AU is an artifact of our choice of a sharp outer boundary for the Oort cloud 
at this point (§4.3.1). All but 2% of the new comets have original energies in the range < x < 10~ 4 AU -1 , 



consistent with our assumption that observed comets in this range are new comets. The mean energy 
of the V\ comets is (l/a) — 3.3 ± 1 x 10~ 5 AU -1 , in good agreement with Heisler's (1990) estimate of 
3.55 x 10~ 5 AU -1 outside of showers. Heisler's Monte Carlo simulations included both the Galactic tide 
and passing stars; the agreement confirms that our omission of stellar perturbers does not strongly bias the 
distribution of new comets. 

The curve in Fig. pd| b shows an analytical approximation to the expected flux of new comets when the 
loss cylinder is full (Wiegert 1996). The agreement between the analytical curve and the distribution of V\ 
comets for a > 30 000 AU confirms that the inner edge of the distribution of new comets is caused by the 
emptying of the loss cylinder as the semimajor axis decreases. The source of the smaller peak at 47 000 AU 
is unclear: if the sample is split into two parts, it appears only in one half, and thus may be a statistical 
fluke even though the deviation from the analytical curve is several times the error bars. In any event it is 
unlikely to play a significant role in determining the overall distribution of LP comets for two reasons: first, 
only a few percent of the V\ comets are involved in the peak; and second the subsequent planet-dominated 
evolution of the V\ comets is relatively insensitive to the comets' original semimajor axes. 

The distributions of perihelion and angular orbital elements for the V\ comets are shown in Fig. |i~2] , 
which can be compared to the observed distributions in Fig. ^. The observed perihelion distribution is 
strongly affected by selection effects, so no comparison is practical there. The angular element distributions 
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Figure 7: The sine of the aphelion latitudes of the 681 LP comets in the ecliptic (a) and Galactic (6) 
reference frames. The heavy line indicates a spherically symmetric distribution. Data taken from Marsden 
and Williams (1993). 



End-state 


Ejection 


Exc. Limit 


Large q 


Short Pd. 


Visible 


Total 


Number 


3807 


57 


15023 


32 


1368 


20286 


Fraction 


0.1877 


0.0028 


0.7406 


0.0015 


0.0674 


1.0000 


Minimum t x 


6.80 


17.2 


7.46 


11.7 


7.14 


6.80 


Median t x 


28.7 


152 


35.2 


29.3 


26.8 


33.3 


Maximum t x 


342 


480 


1182 


72.4 


147 


1182 


Minimum m x 


1 


5000 


1 


6 


1 


1 


Median m x 


8 


5000 


5 


387 


5 


6 


Maximum m x 


4799 


5000 


4872 


3432 


2937 


5000 



Table I: The distribution of end-states of the 20 286 Oort cloud comets with minimum perihelia < 40 AU. 
The minimum, median and maximum lifetimes m x and t x are shown in orbital periods and Myr, respectively. 
No comets suffered collisions with the planets or the Sun, or survived for the lifetime of the Solar System. 
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Figure 8: Distribution of orbital elements for the 109 new comets (1/a < 10 _4 AU _1 ): (a) perihelion 
distance; (b) inclination; (c) longitude of ascending node; (d) argument of perihelion. All angular elements 
are measured in the Galactic frame at the aphelion preceding the first apparition. 
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Figure 9: Equal-area plot of the aphelion directions of the 109 new comets in the Galactic frame. 



Planet 


Jupiter 


Saturn 


Uranus 


Neptune 


Total 


Number of comets 


60 


145 


71 


67 


343 


Number of encounters 


210 


317 


109 


93 


729 


Encounters /comet 


3.5 


2.19 


1.53 


1.39 


2.13 


Collisions 

















Captures 

















Min. distance (R T ) 


0.023 


0.043 


0.074 


0.049 


0.023 


Min. distance (Rp) 


16.0 


38.7 


150 


167 


16.0 


Rj(Rp) 


674 


907 


2030 


3510 




Outer satellite (Rp) 


326 


216 


23 


222 





Table II: Planetary encounter data during the pre-visibility stage for the 20 286 Oort cloud comets with 
minimum perihelia < 40 AU. Encounters for the 57 comets in the Exceeded Time Limit end-state are 
included only up to their 5000*^ orbit. "Captures" are considered to occur when the comet has a planetocen- 
tric eccentricity less than unity at planetocentric pericentre. The radius of the planet's sphere of influence 
i?j (Eq. and the semimajor axis of its outermost satellite are also given, in units of the planetary radius 
R p . 
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2xlO- 5 4xl0- 5 6x10-5 8x10-5 0.0001 2x10" 4xl0 4 6xl0 4 



1/a (1/AU) a (AU) 

(a) (b) 

Figure 11: Distribution of original energies x = 1/a and semimajor axes a for the V\ comets. An additional 
28 comets, 2% of the total, have x > AU" 1 . The curve in (b) is an analytical approximation to the 
expected distribution when the loss cylinder is full, derived in Wiegert (1996). 




Figure 12: Distribution of orbital elements for the V\ comets: (a) perihelion distance; (b) inclination; (c) 
longitude of ascending node; (d) argument of perihelion. All angular elements are measured in the Galactic 
frame when the comet passes 200 AU on its inbound leg. 
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are reasonably consistent between the two figures; in particular the w distributions both show peaks in the 
regions where sin 2uj > 0, reflecting the role of the Galactic tide in creating new comets. 

The aphelion directions of the V\ comets are shown in Fig. 13, which can be compared to the observed 
distribution in Fig. ^. The most striking feature in Fig. [l^ is the concentration towards mid-Galactic 
latitudes, again pointing to the importance of the Galactic tide as a LP comet injector. The real distribution 
of aphelion directions is expected to be much clumpier, due to the injection of comets by passing stars; 
however, the number of new comets in Fig. ^| is too small for any reliable comparisons to be made. 



5.1.1 The longest-lived comets 

Although most comets reach one of the end-states within a few orbits (see Table ||) a small fraction survive 
for much longer times: 57 of the 20 286 initial comets in our simulation triggered the Exceeded Orbit 
Limit flag after 5000 orbits. The population of these comets decays only very slowly and their fate cannot 
be determined without prohibitive expenditures of CPU time. The perihelion distances and semimajor axes 
of these comets on their 5000" 1 orbit are indicated in Fig. [l4|. Also shown is the distance at which they 
cross the ecliptic. Most have nodes and perihelia outside Saturn's orbit, where the energy perturbations are 
relatively small. 



5.2 Post-visibility evolution: the standard model 

We now follow the orbits of the V\ comets forward in time until they reach one of the end-states (obviously, 
the visible end-state is disabled in these simulations). Each time one of these comets makes an apparition 
its orbital elements are added to the set of Voo comets. The Voo comets are to be compared to the observed 
distribution of LP comets. 

The errors in the distribution of elements of the Voo comets are not Poisson, as a single comet may 
contribute hundreds or thousands of apparitions. The errors that we quote and show in the figures are 
determined instead by bootstrap estimation (Efron 1982; Press et al. 1992). 

The "standard model" simulation of post-visibility evolution has no fading, and no perturbers except the 
giant planets and the Galactic tide. 



The distribution of end-states for the standard model is shown in Table III . The Exceeded orbit limit 



end-state (§ L4) is invoked after 10 000 orbits for these simulations, but no comets reach this end-state. The 
mean lifetime is 45.3 orbits, compared to 60 predicted by the gambler's ruin model (Eq. |l5| ). Ejection by 
the giant planets is by far the most common end-state (89% of V\ comets). Most of the remaining comets 
(about 8% of the total) move back out to large perihelion distances. Their median energy when they reach 
this end-state is given by 1/a = 4x 10~ 5 AU^ 1 (a = 25 000 AU); in other words these comets have suffered 
relatively small energy perturbations and remain in the outer Oort cloud. 

The distribution of orbital elements of the Voo comets may be parametrized by the dimensionless ratios 
Xi defined in Eq. ||: the ratio of theoretical parameters \&* for the standard model to the observed parameters 
(Eq. I is 

Xi = -^ = 0.075 ±0.011, X 2 = ^ = 4.4 ± 1.2, X 3 = = 0.61 ± 0.13. (46) 

*1 4-2 *3 

The standard model agrees much better with the predictions of the simple gambler's ruin model (X\ = 0.05, 
A3 = 0.58, see eqs. [l6] and [Tt]) than it does with the observations [Xi = 1). 

The perihelion distribution of the Vx comets in the standard model is shown in Fig. [15[ Although the 
figure represents 52 303 apparitions, the error bars — as determined by bootstrap — remain large, reflecting 
strong contributions from a few long-lived comets: over 45% of the apparitions are due to the 12 comets 
that survive for 1000 or more orbits after their first apparition. This figure can be compared to the observed 
perihelion distribution (Fig. ||), which however reflects the strong selection effects favouring objects near the 
Sun or the Earth. We note that not all perihelion passages made by comets after their first apparition are 
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a (AU) nodal radii (AU) 



Figure 14: For the 57 comets that survived 5000 orbits, (a) their perihelion distance q versus semimajor axis 
a, and (b) the distances of their nodes. In (a), triangles are prograde comets, squares, retrograde. 
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visible: in addition to the 52 303 apparitions made by the Voo comets, there were 9561 perihelion passages 
with q > 3 AU. 

Let the total number of comets with perihelia in the range [q,q + dq] be N(q)dq. A linear least-squares 
fit to Fig. [l5| yields N(q) roughly proportional to 1 + (q/1 AU), similar to Everhart's (1967a) earlier estimate 
of the intrinsic perihelion distribution. The perihelion distribution is not flat, as would be expected if the 
distribution were uniform on the energy hypersurface (Eq. Q). The simulations are noisy enough to be 
consistent with any number of slowly varying functions of perihelion over < q < 3 AU, possibly including 
N(q) oc q 1 / 2 , as proposed by Kresak and Pittich (1978). The estimates of the intrinsic perihelion distribution 
of LP comets published by Everhart, by Kresak and Pittich, and by Shoemaker and Wolfe are indicated on 
Fig. 0. 

The original energy distribution of the comets in the standard model is shown in Fig. |l^, at two 
different magnifications, for all 52 303 apparitions. These figures should be compared with the observations 
shown in Fig. |[ As already indicated by the statistic X\ (Eq. |46|) , the standard model has far too many LP 
comets relative to the number of comets in the Oort spike: the simulation produces 35 visible LP comets for 
each comet in the spike, whereas in the observed sample the ratio is 3:1. This disagreement is at the heart 
of the fading problem: how can the loss of over 90% of the older LP comets be explained? 

These simulations allow us to estimate the contamination of the Oort spike by dynamically older comets. 
There are 1368 V\ comets, of which 1340 have 1/a < 10 -4 AU -1 , but a total of 1475 apparitions are made in 
this energy range in the standard model. Thus roughly 7% of comets in the Oort spike are not dynamically 
new. Of course, this estimate neglects fading, which would further decrease the contamination of the Oort 
spike by older comets. 

Figure 17 shows the inclination distribution of the comets in the standard model. There is a noticeable 
excess of comets in ecliptic retrograde orbits: the fraction on prograde orbits is 15875/52303 » 0.3. This 
is inconsistent with observations, which show an isotropic distribution (Fig. ^a), but consistent with the 
predictions of the gambler's ruin model (Eq. 0). 

Figure [l8] shows the distribution of the longitude of the ascending node and the argument of perihelion, in 
the ecliptic frame. The large error bars suggest that the structure in these figures is probably not statistically 
significant. 

The principal conclusion from this analysis is that the standard model provides a poor fit to the observed 
distribution of LP comets. The standard model agrees much better with models based on a on e-di mensional 
rand om walk , suggesting that the basic assumptions of the analytic random- walk models in §3^3 are valid. 
In §§5^-5^ we shall explore whether variants of the standard model can provide a better match to the 
observations. 



5.2.1 Short-period comets from the Oort cloud 

During our simulations only 68 Oort cloud comets eventually become short-period comets, 36 of them after 
having made one or more apparitions as LP comets. The distributions of inverse semimajor axis, perihelion 
distance and inclination for these comets are shown in Fig. In no case is an Oort cloud comet converted 
to a short-period comet in a single perihelion passage: the largest orbit at the previous aphelion has a 
semimajor axis of only 1850 AU. There is a distinct concentration of orbits near zero ecliptic inclination, as 
expected from studies of captures by Jupiter (Everhart 1972), but the concentration is much less than that 
of short-period comets in our Solar System. The prograde fraction is 44/68 ~ 0.65. 

Our simulation corresponds to approximately 450 years of real time (Eq. ^4|). Thus we deduce that 
68/450 ~ 0.15 short-period comets per year arrive (indirectly) from the Oort cloud (in the absence of 
fading). For comparison, on average five new short-period comets are discovered each year (Festou et al. 
1993a); we conclude that the Oort cloud contributes less than 3% of the population of short-period comets, 
and another source, such as the Kuiper belt, is required. Only about 10% of the known SP comet apparitions 
are Halley-family, and thus the Oort cloud may contribute a significant fraction of these objects, though the 
picture is clouded by the multiple apparitions by individual comets in this sample. 
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End-state 


Ejection 


Large q 


Short pd. 


Total 


Number 


1223 


109 


36 


1368 


Fraction 


0.894 


0.080 


0.026 


1.000 


Minimum t x 


0.296 


2.61 


0.014 


0.014 


Median t x 


1.33 


4.62 


0.67 


1.40 


Maximum t x 


31.7 


71.0 


7.94 


71.0 


Minimum m x 


1 


1 


13 


1 


Median m x 


1 


2 


330 


1 


Maximum m x 


5832 


2158 


4277 


5832 



Table III: The distribution of end-states of the V\ comets in the standard model. The minimum, median 
and maximum lifetimes t x of these comets are measured from their first apparition in Myr. No comets suffer 
collisions with the planets or Sun, or survive for the age of the Solar System. 




c (AL) 

Figure 15: Distribution of perihelion distances q for the comets in the standard model. Error bars 
are determined from bootstrap estimators and represent one standard deviation. The curves are Everhart's 
(1967a, dotted line), Kresak and Pittich's (1982, solid line) and Shoemaker and Wolfe's (1982, dashed line) 
estimates of the intrinsic perihelion distribution. The correct normalizations arc unclear, and have been 
made somewhat arbitrarily. 
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Figure 17: Distribution of the cosine of the inclination for the Voo comets in the standard model, (a) at 
perihelion in the ecliptic frame, and (b) at 200 AU on the inbound leg in the Galactic frame. The heavy 
line indicates a uniform distribution. 
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5.2.2 Planetary encounter rates 



Close encounters of the Voo comets with the giant planets are described in Table [V . Note the high frequency 
of multiple encounters between a giant planet and a single comet, though this does not indicate capture by 
the planet in the traditional sense (i.e. planetocentric eccentricity less than unity). 

Since our simulation corresponds to roughly 450 years of real time (Eq. |44|) we can calculate the rate of 
close encounters between the LP comets and the giant planets. During the combined pre- and post-visibility 
phases of the comets' evolution, a total of 253 encounters were recorded for Jupiter, 333 for Saturn, 111 for 
Uranus and 96 for Neptune. These numbers translate to total currents J p of 0.56, 0.74, 0.25 and 0.21 comets 
per year passing through the spheres of influence (Eq. pi| ) of Jupiter through Neptune respectively. 

If we assume that these currents J p reflect a uniform flux of LP comets across the sphere of influence of 
each planet, then the rate of impacts between LP comets and the giant planets can be deduced to be 

M P Y 4/5 {Rp\ 2 (. 2 M, 



where a p and R p are the planets' semimajor axis and radius, M p is the planetary mass, and the second 
term is a crude correction for gravitational focusing, assuming the comets are on nearly parabolic orbits. 
The resulting collision rates are 1.0 x 10~ 5 , 5.0 x 10~ 6 , 2.4 x 10~ 7 , 1.3 x 10 -7 per year for Jupiter through 
Neptune respectively. It should be noted that Comet Shoemaker-Levy 9, which collided with Jupiter in July 
of 1994, was not a LP comet but rather a Jupiter- family comet (Benner and McKinnon 1995). 

5.3 Post-visibility evolution: the effect of non-gravitational forces 

Asymmetric sublimation of volatiles leads to significant non-gravitational (NG) forces on comets. As de- 



scribed in § |4.2] , we specify NG forces using two parameters A\ and A 2 . The parameter A\ is proportional to 
the strength of the radial NG force, and is always positive, as outgassing accelerates the comet away from the 
Sun. The parameter A2 is proportional to the strength of the tangential force, is generally less than A\, and 
may have either sign depending on the comet's rotation. Comet nuclei are likely to have randomly oriented 
axes of rotation, with a corresponding random value of A2. Rather than make a complete exploration of the 
available parameter space for A\ and A2, we shall investigate a few representative cases. 
We assume that \A 2 \ = O.lAi, and consider two distributions for the sign of A 2 : 

1. Half the comets have positive values of A 2 , half negative, and the sign of A 2 is constant throughout 
a comet's lifetime — as if the axis of rotation of the nucleus remained steady throughout the comet's 
dynamical lifetime. 

2. The sign of A 2 is chosen at random after each perihelion passage — as if the axis of rotation changed 
rapidly and chaotically. 

We examined four values of A\: 10~ 8 , 10~ 7 , 10~ 6 , 10 -5 AU day -2 . The first two of these are reasonably 
consistent with the NG forces observed in LP comets (Marsden et al. 1973). The two remaining values for 
Ai are probably unrealistically large. 

Figure |2^ and Table ^ illustrate the effects of NG forces on the energy and perihelion distributions, 
and on the parameters Xi defined in Eq. ^, which should be unity if the simulated and observed element 
distributions agree. The figure shows that NG forces do decrease the number of dynamically older comets 
relative to the number of new comets and hence improve agreement with the observations (i.e. increasing 
Ai, decreasing X 2 ); however, the same forces erode the population of comets at small perihelion distances, 
thereby worsening the agreement with the observed perihelion distribution. Even unrealistically large NG 
forces cannot bring the distribution of inverse semimajor axes into line with observations, and these produce 
an extremely unrealistic depletion of comets at small perihelia. 

The effects of NG forces can be summarized as follows: 
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Planet 


Sun 


Jupiter 


Saturn 


Uranus 


Neptune 


Total 


Number of comets 


7 


28 


12 


2 


3 


52 


Number of encounters 


16 


43 


16 


4 


3 


82 


Encounters / comet 


2.3 


1.5 


1.3 


2.0 


1.0 


1.6 


Collisions 




















Captures 



















Min. distance 




0.018 


0.086 


0.17 


0.16 


0.018 


Min. distance (R p ) 


1.61 


12.5 


77.9 


335 


553 


12.5 


Outer satellite (R p ) 




326 


216 


23 


222 





Table IV: Planetary and solar encounter data for the Voo comets in the standard model (post-visibility stage) . 
The distance to each planet's outermost satellite is given in the last row. 
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Figure 19: The distribution of the inverse semimajor axis 1/a, perihelion distance q and cosine of the ecliptic 
inclination i for the short-period comets originating in the Oort cloud. The distribution of 1/a on the left is 
measured at the aphelion previous to, and the other distributions measured at, the initial perihelion passage 
as a short-period comet. 
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The semimajor axis perturbation due to radial NG forces averages to zero over a full orbit (assuming 



that the radial force is symmetric about perihelion, as in the model discussed in £4.2). Thus radial 
forces have little or no long-term effect on the orbital distribution. 

• Positive values of the tangential acceleration A 2 reduce the tail of the population, resulting in an 
increase in X\ towards unity and improving the match with observations, but erode the population at 
small perihelia, a depletion which is not seen in the observed sample. 

• Negative values of A 2 preserve a reasonable perihelion distribution, but increase the number of comets 
in the tail of the energy distribution, thus reducing X\ so that the disagreement between the observed 
and simulated energy distribution becomes even worse. 

We have also conducted simulations with a more realistic model for observational selection effects (Eq. |]) 
but this does not alter our conclusions. 

Although we have not exhaustively explored the effects of NG forces on the LP comet distribution, we 
are confident that conventional models of NG forces cannot by themselves resolve the discrepancy between 
the observed and predicted LP comet distribution. 



5.4 Post-visibility evolution: the effect of a solar companion or disk 

In this section we investigate the influence of two hypothetical components of the Solar System on the 
evolution of LP comets: 

1. A massive circumsolar disk extending to hundreds of AU or even further. Such a disk might be an 
extension of the Kuiper belt (Jewitt et al. 1996) or related to the gas and dust disks that have been 
detected around stars (especially 3 Pictoris) and young stellar objects (Ferlet and Vidal-Madjar 1994). 
Residuals in fits to the orbit of Hallcy's comet imply that the maximum allowed mass for a disk of 
radius r is roughly (Hogg et al. 1991) 

M max .10M e (^) 3 . (48) 

Current estimates of the mass in the Kuiper belt are much smaller, typically ~ O.IM© from direct 
detection of 100 km objects (Jewitt et al. 1996) or from models of diffuse infrared emission (Backman 
et al. 1995), but these are based on the uncertain assumption that most of the belt mass is in the 
range 30-50 AU. The disk around (3 Pic is detected in the infrared to radii exceeding 1000 AU (Smith 
and Terrile 1987); the dust mass is probably less than 1M© (Artymowicz 1994), but there may be 
more mass in condensed objects. 

2. A solar companion, perhaps a massive planet or brown dwarf, orbiting at hundreds of AU. Residuals 
in fits to the orbits of the outer planets imply that the maximum allowed mass for a companion at 
radius r is roughly (Tremaine 1990; Hogg et al. 1991) 

m -^ 100M -(mau) 3 - < 49 > 

There are also significant but model-dependent constraints on the characteristics of a solar companion 
from the IRAS infrared all-sky survey (Hogg et al. 1991). 

To reduce computational costs, we used the V\ comets as a starting point for these investigations; that is, 
the effect of the disk or companion is ignored before the comet's first apparition (more precisely, we started 
the integration at the aphelion preceding the comet's initial apparition, in order to correctly calculate any 
perturbations occurring on the inbound leg). Starting at this point is an undesirable oversimplification, but 
one that should not compromise our conclusions. 
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Figure 20: Distribution of the inverse semimajor axis l/o and perihelion distance q for Vac comets subjected to 
non-gravitational forces. Left panels: constant values for A2, half positive, half negative. Right panels: the sign of 
A2 is randomised for each perihelion passage. From the top down, Ai = 1(T 8 , 10 -7 , 1(T 6 and 10~ 5 AU day" 2 , 
with \A2 1 = O.lAi. The bottom line of panels is for comparison, and includes the standard model (left side) and 
the observations (right side). The observed perihelion distribution includes curves indicating the estimated intrinsic 
distribution (see Fig. for details) . 
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5.4.1 Circumsolar disk 

The circumsolar disk is represented by a Miyamoto-Nagai potential (e.g. Binney and Tremaine 1987), 

V disk (x, y, z) = : — 7T7T/2 - ( 50 ) 



+ y 2 + (a d + y^z 2 + b 2 d y 



Here Md is the disk mass, and a d and bd are parameters describing the disk's characteristic radius and 
thickness. We assume that the disk is centered on the Solar System barycenter and coplanar with the 
ecliptic. We considered disk masses Md of 0.1, 1, and 10 Jupiter masses, disk radii ad of 100 and 1000 AU, 
and a fixed axis ratio bd/a d = 0.1. The two more massive disks with ad — 100 AU are unrealistic because 
they strongly violate the constraint (^8|) , but we examine their effects in order to explore as wide a range of 
parameters as possible. 

Comets arriving from the Oort cloud have fallen through the disk potential and hence are subjected to a 
shift in their original inverse semimajor axis. This offset can be as large as 2 x 10~ 4 AU -1 for a 10 Jupiter- 
mass disk with radius 100 AU, but is much smaller for disks that do not already violate the observational 
constraint ([l8|). This shift is not shown in the Figures below, for which the semimajor axis is measured at 
aphelion. 



As usual, the Perihelion Too Large end-state (§ [L4|) was entered if q > 40 AU and sin2o) > 0. The 
assumption that such comets are unlikely to become visible in the future is only correct if the torque is 
dominated by the Galactic tide, and this may not be the case when a disk is present. However, there is no 
significant difference in the numbers or semimajor axes of the comets reaching this end-state in simulations 
with and without a circumsolar disk, suggesting that evolution to this end-state is indeed dominated by the 
Galaxy. 

The results from simulations including a circumsolar disk are displayed in Fig. |l| and Table [v| One 



plot of the energy distribution in Fig. 21 shows a strong peak near 1/a = 0.02 AU -1 ; as the large error bars 
suggest, this peak is caused by a single comet and has little statistical significance. 

The principal effect of the disk is to exert an additional torque on the comets, resulting in oscillations of 
the comet's perihelion distance. This effect normally increased the comet's lifetime, as the risk of ejection 
is greatly reduced when the comet is outside Saturn's orbit. The perihelion oscillations also enhance the 



probability of collision with the Sun (Tabic VI) 



The perihelion distribution of visible comets is not strongly affected by the disk. The presence of a massive 
disk reduces the number of dynamically old comets (because their perihelia are no longer nearly constant, 
only a fraction of them are visible at any given time), but not enough so that the energy distribution is 
consistent with the observations. This conclusion is confirmed by examining the X parameters in Table |V]| , 
which should be unity if the simulated element distribution agrees with the observations (cf. Eq. ||). The 
values of X\, which measures the ratio of number of comets in the spike to the total number, are far smaller 
than unity even for the most massive disks. Increasing the disk mass tends to improve X2 and X3 for the 
1000 AU disk, but degrades the fit for the 100 AU disk. There is no set of disk parameters than comes close 
to producing a match with observations. Using a more elaborate model for selection effects (Eq. [[]) does not 



alter this conclusion (see bottom half of Table VI ) 



We conclude that a circumsolar disk cannot by itself resolve the discrepancy between the observed and 
predicted LP comet distribution. 



5.4.2 Solar companion 

For simplicity, we shall assume that the solar companion has a circular orbit in the ecliptic (the orientation 
and eccentricity of the companion orbit should not strongly affect its influence on the LP comets since the 
comets are on isotropic, highly eccentric robits). 

We examined companion masses Mx of 0.1, 1 and 10 Jupiter masses and orbital radii of 100 and 1000 AU. 
The most massive companion at 100 AU is unrealistic because it strongly violates the constraint (|^). As in 
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A 2 


Total 
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Tail 


Prograde 


Xi 


x 2 


X 3 


(m) 


0.0 


0.0 


52303 


1473 


15004 


15875 


0.07 


4.37 


0.61 


45.4 


1.0 


0.1 


35370 


1457 


7368 


12381 


0.11 


3.17 


0.70 
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1.0 


-0.1 


57819 
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19364 
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0.07 
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0.73 


51.0 


1.0 


±0.1 a 


44383 


1461 


13705 


19021 


0.09 


4.70 


0.86 


38.4 
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45899 


1425 


16628 


18504 


0.08 


4.42 


0.80 
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±10 a 


30660 


1341 


11296 


11012 


0.12 


5.61 


0.72 


33.1 
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±100 a 


13248 
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5432 
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6.24 


0.88 


14.4 


1.0 


±0.1 b 


49642 


1450 


13203 


16387 


0.08 


4.05 


0.66 


46.7 


10 


±l b 


45202 


1448 


13631 


17311 


0.08 


4.59 


0.76 


41.4 


100 


±10 h 


25774 


1364 


4969 


11452 


0.14 


2.93 


0.88 


27.7 


1000 


±100 b 


9878 


1035 


1536 


5042 


0.28 


2.37 


1.02 


13.2 



Table V: Parameters of the distribution of Vqo comets subjected to non-gravitational forces. The superscript 
a indicates that half the sample have positive A 2 , half negative; b indicates that A 2 has a randomly chosen 
sign for each perihelion passage. "Total" is the total number of apparitions (i.e. perihelion passages with 
q < 3 AU), "Spike" is the number of these with original inverse semimajor axes 1/a < 10~ 4 AU -1 , "Tail" is 
the number with 0.0145 AU -1 < 1/a < 0.029 AU -1 , and "Prograde" the number with ecliptic inclination less 
than 90°. The parameters Xi are defined in Eq. [| The mean lifetime in orbits (m) includes all perihelion 
passages, whether visible or not, after the initial apparition. The units of A\ and A 2 are AU day -2 . 
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18769 
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Table VI: Parameters of the distribution of V^, comets, when the Solar System contains a circumsolar disk. 
The disk mass Md is measured in Jupiter masses and the disk radius ad is measured in AU. The rightmost 
column indicates the number of comets that collided with the Sun. The superscript d indicates that the 
discovery probability from Eq. [I] has been applied. The definitions of the other columns are the same as in 
Table [v| 



5 RESULTS 



40 





I/O (1/AU) " q (AU) |/a (1/AU) q (AU) 



Figure 21: Distribution of the inverse scmimajor axis 1/a and perihelion distance q for the Vqo comets, when 
the Solar System contains a massive circumsolar disk. Left panels: characteristic disk radius ad = 100 AU. 
Right panels: disk radius ad = 1000 AU. From the top down, the disk masses are 0.1, 1 and 10 Jupiter 
masses. The bottom line of panels is for comparison, and includes the standard model (left side) and the 
observations (right side). The observed perihelion distribution includes curves indicating the estimated 
intrinsic distribution (see Fig. fHi for details). 
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the previous subsection, the original semimajor axes of the comets are measured at aphelion, and thus do 
not include the energy offset caused by their fall through the companion's gravitational potential. 



The results are presented in Fig. 22 and Tabic VII. The X parameters are listed in Table VII. As the 
companion mass is increased, the fraction of prograde to total comets (A3) improves. However, X\ and Xi 
remain far from unity. There is no evidence that a solar companion can significantly improve the agreement 
between the observed and predicted LP comet distribution. 

5.5 Post- visibility evolution: fading 



The concept of fading was introduced in §§3.3 and 3.4. We use the term "fading" to denote any change in 



the intrinsic properties of the comet that would cause it to disappear from the observed sample. Our focus 
is on modeling the fading process empirically, rather than attempting to elucidate the physical processes 
involved. The distributions of inverse semimajor axis and ecliptic inclination will serve as our primary 
fading benchmarks, through the values of the parameters Xi, X 2 and X 3 (Eq. ||). 

We shall generally assume that fading depends only on the number of apparitions (perihelion passages 
with q < 3AU). We parametrize the fading process by a function $ m (cf. Eq. ^J), the probability that a 
visible new comet survives fading for at least to apparitions (thus $1 = 1). 



We shall conduct simulations with and without plausible non-gravitational (NG) forces (§§4.2, 5.3). When 
NG forces are included, we shall use the parameters A\ = 10~ 7 AU day -2 , A2 = ±10~ 8 AU day" 2 , A3 = 0, 
with a random sign for A2 at each perihelion passage (henceforth the "standard NG model" ) . 

The most direct way to determine the fading function <E> m would be to break down the simulated data set 
into individual distributions, one for each perihelion passage i.e. {V\, V2, V3, . . .}, and then fit the observed 
distribution of orbital elements to the parameters $1, $2, ■ ■ • where $ m +i < $ m - Unfortunately, this problem 
is poorly conditioned. Instead, we shall experiment with a few simple parametrized fading functions. 

5.5.1 One-parameter fading functions 

The fading functions we shall examine include: 

a) Constant lifetime Each comet is assigned a fixed lifetime, measured in apparitions. Thus 

$ m = 1, to < m v , $ m = 0, to > m v . (51) 

b) Constant fading probability Comets are assigned a fixed probability A of fading, per apparition. Thus 

$ ro = (l-A) m - 1 . (52) 

c) Power-law The fraction of comets remaining is 

$ m = m- K , (53) 

where k is a positive constant. 

We have also investigated fading functions in which $ depends on the elapsed time t since the first 
apparition. Such laws are less physically plausible than fading functions based on the number of apparitions, 
since by far the harshest environment for comets occurs as they pass perihelion; and in fact the functions 
<!>(£) that we investigated all produced relatively poor matches to the observations. Also, fading functions 
in which $ depends on the number of perihelion passages produce results very similar to laws based on the 
number of apparitions. 

The results from the fading laws (|5l|)-(p3|) are shown in Figs. 23 to [2^. The first of these figures displays 



the X parameters assuming LP comets have a constant lifetime in apparitions (model [a]). The presence or 
absence of NG forces (bottom vs. top panels), or the use of two different visibility criteria (left vs. right 
panels) has very little effect on the results. The spike/total ratio matches observations (i.e. X\ — 1) at 
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l/a (1/AU) ' q (AU) 1/0 (1/AU) q (AU) 

Figure 22: Distribution of the inverse scmimajor axis l/a and perihelion distance q for the Vqo comets, when 
the Solar System contains a massive solar companion. Left panels: companion orbital radius of 100 AU. 
Right panel: orbital radius 1000 AU. From the top down, the companion masses are 0.1, 1 and 10 Jupiter 
masses. The bottom line of panels is for comparison, and includes the standard model (left side) and the 
observations (right side). The observed perihelion distribution includes curves indicating the estimated 
intrinsic distribution (see Fig. fHi for details). 
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Figure 23: The values of the parameters A, for a fading function with a fixed lifetime of m v apparitions 
(model [a], Eq. If the simulation agrees with the observations then Xj = 1, i = 1,2,3. The parameter 
X\ is based on the fraction of LP comets in the Oort spike (solid curve); Xi is based on the fraction of 
comets in the energ y ta il, x > 0.0145 AU -1 (dotted curve); X3 is based on the fraction of prograde comets 



(dashed curve) (cf. §2.8). The panels on the left are based on the visibility criterion q < 3 AU, and those on 
the right are based on the visibility probability (Eq. |l|). The upper panels are based on the standard model 
with no NG forces, and the lower panels are based on the standard NG model. 
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Figure 24: The values of Xi given a fixed fading probability A per apparition (model [b], Eq. j52"l) . 
details see the caption to Fig. 
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m v ~ 10, but the tail/total ratio is far too low at that point (X 2 <C 1)- The tail/total ratio is right at 
m v ~ 100, but Xi is now too low. The ratio X3 is typically close to but below unity. The model does not 
match the observations for any value of the parameter m v . 

Figure [24] displays the behaviour of the parameters Xj given a fixed fading probability A per apparition 
(model [b]). Once again, the results are almost independent of NG forces and the visibility criterion, and 
there is no value for the parameter A that matches the observations (X; = 1). 

Figure |2^ shows the parameters Xj for a power-law fading function (model [c]). Although the match 
is not perfect, an exponent k — 0.6 ± 0.1 provides a much better match than the previous two models: 
Xi = 0.73 ± 0.09, X 2 = 0.96 ± 0.26, and X 3 = 0.95 ± 0.12 when the standard NG model and discovery 
probability (Eq. [l]) are used. The distributions of orbital elements are shown in Fig. ^6|, to be compared 
with the observed distributions in § ^. For m ^> 1 this fading law is the same as an empirical law suggested 
by Whipple (1962), cf> m = $ m — $ m+ i oc m^ K_1 ; Whipple estimated k — 0.7. 



5.5.2 Other fading functions 

We have also examined several two-parameter fading functions: 

d) Two populations Suppose that the Oort cloud contains two populations of comets, distinguished by 
their internal strength. The first and more fragile set is disrupted after m v apparitions, while the more 
robust comets, comprising a fraction / of the total, do not fade at all. Thus 

$ m = 1, m < m v , $ m = /, m > m v . (54) 



e) Constant fading probability plus survivors One population has a fixed fading probability A per 
apparition, while the more robust comets, comprising a fraction / of the total, do not fade at all. Thus 

$ m = (l-/)(l-Ar- 1 + /. (55) 



f) Offset power law The fading function is chosen to be 

<& m = [(m + /?)/(l + /?)]-*, (56) 

The results of model (d) are shown in Fig. |27]. In most cases the fit is worse than in the one-parameter 
model (a), shown by the heavy lines, because the prograde fraction described by X3 is lower when some 
of the comets do not fade. The best match is for the standard NG model with visibility probability (|l|) 
(lower right panel). Here the parameters m v = 6, / = 0.04 yield Xi = 0.83 ± 0.10, X 2 = 0.91 ± 0.26, 
X3 = 0.96 ± 0.11, slightly better than the match for model (a). This model is reminiscent of Weissman's 
(1978) favoured model, in which 85% of LP comets had significant fading probabilities while the reminder 
survived indefinitely. 

Model (e) is a generalization of the one-parameter model (b) but ordinarily does no better: the match to 
observations is usually best when the survivor fraction / is set to zero, and gets worse as / increases. Model 
(f) also does no better than its one-parameter counterpart, model (c). 

Finally, we examine 



g) Other published fading functions: In § 3.3, we described a number of fading functions deduced in 
previous studies. Oort (1950) took ip x = 0.8, ip m = 0.014 for m > 1; Kendall (1961) took ipi = 0.8, 
ip m = 0.04 for m > 1; Whipple (1962) took 4> m oc m" 17 ; Weissman (1978) took / = 0.15, A = 0.1 (cf. 
Eq. |55|); Everhart (1979) took $1 = 1, = 0.2 for m > 1; Bailey's (1984) fading law is described by 
Eq. (|20|); and Emel'yanenko and Bailey (1996) assume <f> m = 0.3 but add a probability k* = 0.0005 



that the comet is "rejuvenated" . In Tabic Vlllj , we have listed the values of Xj obtained for all these 



fading models (the results in the Table are based on the model that includes the discovery probability 
([l]) and standard NG forces; other models give very similar results). Many provide reasonable matches 
to the data but none do as well as our best fits. 
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M x 


ax 


Total 
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Tail 


Prograde 


Xi 


x 2 
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Rq 


0.1 


100 


40662 


1451 
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14074 


0.11 


3.07 


0.67 


43.1 


1 


0.1 
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10057 
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0.09 
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0.53 
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1 
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38397 
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0.12 
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0.47 
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4 


1 
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35940 
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13544 


0.12 


3.56 


0.73 
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1 


10 
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1379 


3365 


5846 


0.28 


3.10 


0.76 


66.0 
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10 


1000 


28600 


1400 
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0.15 


3.92 


1.05 


146.3 


2 


0.1 d 
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25300 


944 


6762 


8893 


0.11 


3.66 


0.68 


43.1 


1 


0.1 d 


1000 


31376 


975 


6206 


8623 


0.09 


2.71 


0.53 


44.4 


1 


l d 


100 


27918 


963 


4764 


6047 


0.10 


2.34 


0.42 


85.4 


4 


l d 


1000 


24740 


943 


6713 


8281 


0.12 


3.72 


0.65 


68.1 


1 


10 d 


100 


9749 


928 


2197 


4059 


0.29 


3.09 


0.81 


66.0 


4 


10 d 


1000 


22177 


1030 


6052 


12649 


0.14 


3.74 


1.11 


146.3 
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Table VII: Parameters of the distribution of Voc comets, when the Solar System contains a massive solar 
companion. The companion mass Mx is in Jupiter masses, and its orbital radius ax is measured in AU. The 
rightmost column indicates the number of comets that collided with the Sun. The superscript d indicates 
that the discovery probability from Eq. [j] has been applied. The definitions of the other columns are the 
same as m Table |V} 



Name 


Xi 


x 2 


X 3 


Oort 


0.66 ±0.09 


1.21 ±0.44 


0.92 ±0.13 


Kendall 


0.99 ±0.12 


0.59 ±0.23 


0.92 ±0.11 


Whipple 


0.97 ±0.11 


0.58 ±0.16 


0.95 ±0.11 


Weissman 


0.50 ±0.07 


2.07 ±0.58 


0.97 ±0.14 


Everhart 


0.47 ±0.07 


2.60 ±0.72 


0.97 ±0.15 


Bailey 


0.82 ±0.11 


1.68 ±0.63 


1.07 ±0.13 


Emcl'yanenko 


0.69 ±0.08 


0.16 ±0.05 


0.94 ±0.10 



Table VIII: The values of Xi for the preferred fading models of Oort (1950), Kendall (1961), Whipple 
(1962), Weissman (1978), Everhart (1979), Bailey (1984), and Emel'yanenko and Bailey (1996). The results 
are based on the model which includes the discovery probability (pi) and the standard NG forces. 
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Figure 26: The distribution of the inverse semimajor axis 1/a, perihelion distance q and cosine of the ecliptic 
inclination i for a power-law fading function with exponent k = —0.6 (Eq. |53| ). These simulations are based 
on the standard NG model and the visibility probability in Eq. [l[ 
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Figure 27: The values of Xi given a two-parameter fading function in which a fraction 1 — / survives for m v 
apparitions, while a fraction / survives forever (model [d], Eq. (EJ). The fractions / for the different curves 
are 0, 0.01, 0.04, 0.07, 0.1, 0.2, 0.4, 0.6, 0.8 and 1, beginning with the heavy lines. For further details see 
the caption to Fig. |2^. 
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6 Summary 

The LP comets provide our only probe of the properties of the Oort comet cloud. The expected distribution 
of their orbital elements is only weakly dependent on the properties of the Oort cloud and is straightforward — 
though not easy — to predict if the distribution is in a steady state. Thus a central problem in the study of 
comets is to compare the predicted and observed distributions of the orbital elements of the LP comets. 

We have simulated the dynamical evolution of LP comets from their origin in the Oort cloud until the 
comets are lost or destroyed. We have integrated the comet trajectories under the influence of the Sun, the 
giant planets, and the Galactic tide. In some cases we have included the effects of non-gravitational forces, 
a hypothetical circumsolar disk or solar companion, and the disruption or fading of the comet nucleus. We 
have not included the effects of passing stars on the Oort cloud; these add a random component to the 
expected distribution of LP comets which is more difficult to model but is not expected to strongly affect the 



distribution except during rare comet showers (cf. §4.1.3). Our conclusions from these simulations include 
the following: 

The Oort cloud presently contains roughly 5 x 10 11 ($ nGW /3 yr -1 ) objects orbiting between 10 000 and 
50 000 AU from the Sun (Eq. ||), assuming that the cloud is in a steady state and that the number density 
in the cloud is proportional to r~ 35 (Duncan et al. 1987); here $ n ew is the observed current of new comets 
with perihelion < 3AU. This estimate depends strongly on uncertain assumptions about the density and 
extent of the inner Oort cloud; a more reliable parameter is that the number of comets in the outer Oort 
cloud (a > 20 000AU) is 4 x 10 n ($ncw/3 yr" 1 ). 



Over 90% of the comets in the Oort spike (1/a < 10 4 AU are making their first apparition (§5.2), and 



only 2% of new comets have energies outside the spike (§5.1). The Oort cloud provides only a few percent 
of the observed short-period comets, and even fewer if LP comets fade. Thus another source, such as the 
Kuiper belt, must provide the bulk of the short-period comets. On the other hand, a significant fraction of 
the Halley-family comets may arise in the Oort cloud; however, biases in and the small size of our both the 
observed and simulated Halley-family comets render this estimate very approximate. 

LP comets collide with Jupiter and Saturn roug hly once per 10 5 yr if $ ncw = 3 vr" 1 (§ |5~2^ ). 

This research does not explain the existence of comets on hyperbolic original orbits (see Fig. |l|). The ex- 
cess velocities are small, corresponding to roughly — 10~ 4 AU -1 in inverse semimajor axis, but are larger than 
those produced by the Galactic tide (~ — 10~ 6 AU" 1 ), by plausible non-gravitational forces (~ — 10~ 5 AU -1 ) 
or by a circumsolar disk or solar companion small enough to be compartible with the distribution of bound 
orbits. 



Using simple models based on a one-dimensional random walk (§3.3), many investigators, starting with 
Oort (1950), have concluded that the observed energy distribution of LP comets is incompatible with the 
expected steady-state distribution, unless most new comets are destroyed before their second or subsequent 
perihelion passage. We have shown that this "fading" problem persists in a simulation that follows the comet 
orbits in detail. 

Non-gravitational forces play a significant role in shaping the distributions of the orbital elements of the 
LP comets, but are too small by at least two orders of magnitude to resolve the fading problem (§ |4.2| ). Hy- 
pothetical additional components of the Solar System such as a massive circumsolar disk or solar companion 



also do not resolve the fading problem (§5.4) 



We can match the observed distribution of orbital elements to the expected steady-state distribution with 
at least two fading functions: (a) a one-parameter power-law (Eq. ^3|) with exponent k ~ 0.6 (Whipple 1962); 
(b) a two-population model (Eq. p4| ) in which approximately 95% of comets survive for roughly six orbits and 
the remainder do not fade (the latter model is also roughly consistent with the observed splitting probabilities 
of dynamically new LP comets, approximately 0.1 per orbit; see Weissman 1980). The observation that the 
cratering rate is roughly compatible with the rate expected from the current known populations of comets 



and asteroids (§ 3.4) suggests — within the large uncertainties — that fading occurs through the fragmentation 
or disruption of the comet nucleus rather than through the production of a single "dead" body. We also note 
the lack of strong fading in new comets during single perihelion passage, which might be expected if fading 
were due to loss of volatiles from an intact nucleus. One possible model compatible with these observations 
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is that the nucleus of a new comet fragments during its first apparition, but that the separation velocity of 
the fragments is small enough that they remain within the coma until well past observability. 

Although physically plausible, fading remains an ad hoc explanation for the distribution of LP comet 
orbits which has not been independently confirmed, and we should remain alert for other possible explana- 
tions. 

We thank Tom Bolton, Ray Carlberg, Martin Duncan, Bob Garrison and Kim Innanen for helpful 
discussions and advice. This research was performed at the University of Toronto and the Canadian Institute 
for Theoretical Astrophysics, and has been supported in part by the Natural Sciences and Engineering 
Research Council of Canada. 
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